Link Copied
A quick and dirty attempt at the Student's t Test for million element sets, where I pick out the first peak between indices of the count array from 89 to 101 provides a set test data to determine if the natural frequency changes with time. It appears to - each million elements is 100 days and each million elements contains about 4% density in the required range.
subroutine mean(i,j,countFZ, meanA, stdevA,countFZA,meanFZA, standevFZA, flag)
implicit none
integer i, j, k, l
integer countFZ(kl)
integer countFZA(kl, mn_9)
REAL (KIND=dp) meanFZA(mn_9),standevFZA(mn_9)
integer flag
REAL (KIND=dp) counter, count, counterV, countold,t, sp,sp1
REAL (KIND=dp) meanA, variance, stdevA, meanold, varianceold
counter = 0.0
count = 0.0
counterV = 0.0
meanA = 0.0
variance = 0.0
do 100 k = i,j
counter = counter + real(countFZ(k)) * real(k)
count = count + real(countFZ(k))
! write(*,1000)k, countFZ(k), counter, count
1000 Format(i7,3(' ',f12.0) )
100 end do
meanA = real(counter)/real(count)
count = 0.0
do 101 k = i,j
counterV = counterV + (real(countFZ(k))) * ((real(k) - meanA) ** 2.0)
count = count + real(countFZ(k))
! write(*,1001)k, countFZ(k), counter, counterV, count
1001 Format(i7,4(' ',f12.0) )
101 end do
variance = real(counterV)/real(count)
write(*,500)0, meanA, variance, sqrt(variance), count
meanold = meanA
varianceold = variance
countold = count
500 Format(" Data set :: ", i4,3(' ',f10.3), ' Count :: ',f7.0, ' Sp : ', f10.3, ' sp1 : ',f10.3,' t-Stat : ',f10.3)
stdevA = sqrt(variance)
write(*,120)flag
120 Format(" Flag for the data size :: ", i3)
do 300 l = 1, flag
counter = 0
count = 0
do 200 k = i,j
counter = counter + real(countFZA(k,l)) * real(k)
count = count + real(countFZA(k,l))
! write(*,1000)k, countFZA(k,l), counter, count
200 end do
meanFZA(l) = counter/count
count = 0.0
counterV = 0.0
do 400 k = i,j
counterV = counterV + (real(countFZA(k,l))) * ((real(k) - meanFZA(l)) ** 2.0)
count = count + real(countFZA(k,l))
! write(*,1001)k, countFZA(k,l), counter, counterV, count
400 end do
variance = counterV/count
standevFZA(l) = sqrt(variance)
sp = sqrt((((countold-1.0)*varianceold)+((count-1.0)*variance))/(countold+count-2))
sp1 = sp*sqrt((1.0/countold)+(1.0/count))
t = (meanFZA(l) - meanold)/sp1
write(*,500)l,meanFZA(l),variance,sqrt(variance),count,sp,sp1,t
meanold = meanFZA(l)
varianceold = variance
countold = count
300 end do
return
end subroutine
I took the standard algorithm from wikipedia -- not equal means or variances.
A t-Stat of 31 to 80 represents an almost zero chance they are the same - mean is 95 column and the variance is the 1.15 column. The counts make is so accurate.
Got to check it and put in error checking.
There is a temperature dependence on the stiffness, but it is impossible to discern at 100 day blocks with any accuracy.
Pavel:
Thank you for the explanation.
The data is in a format where is is simpler to take the first moment and then the second moment to get the variance.
The interesting issue is the amplitude of the t-Stat -- it is way beyond the published tables. So I will need to encode estimating the probability as well. A 1:1000 p value has a t-Stat of about 3.5 for this number of elements, but we are at 37. Interestingly the same problem arises in EXCEL with large data sets when you start looking for 5 and 6 standard deviations p estimates.
I am pretty sure you have the Box Mueller Random Number generator. There is a standard reference that some have used for discussing the Box Mueller Theory from the 1950s. When I was writing a paper, I thought I should get a copy of this original paper if I wanted to reference it - you are supposed to read what you reference. It came from a famous statistical group at a major Eastern University. I sent them a message and asked for a pdf of the paper as I could not find it online. I wanted to see the original math, so I could follow the Fortran.
They came back and said they had no idea, but they finally tracked down a copy and scanned it for me. Apparently I was the first person to ask for it in a long time. There is a brief version in a journal in 1958, but not the entire mathematics, if I remember correctly.
Everyone assumes that someone else has seen it so it is ok to reference.
The same thing happens with a famous book on Seismology in Japan. This book is quoted often, a friend once found there is one acknowledged copy of this book in a Scottish University and no others. There are no pdf's.
