Using nonlinear regression
Do you want to estimate the parameters of a lognormal distribution or learn how to do nonlinear regression in R? If the former, as far as I know, the best known method is maximum likelihood, for which the answer is to compute mean and standard deviations of the logs. This assumes you are talking about the 2-parameter lognormal. I don't know the best method for a 3-parameter lognormal. If that's what you want, PLEASE do read the posting guide! "http://www.R-project.org/posting-guide.html". Doing so can increase the chances of getting a useful reply. If you want examples of nonlinear regression, have you considered "nls" and "optim"? spencer graves
Mark Miller wrote:
Attached is a copy of my code, the data and the plots obtained by varying terms manually, I was told that nonlinear regression in R could find the values for me, but I am unable to figure how exactly I could implement this. Any help would be very greatly appreciated as I am completely stuck on this problem. On Thursday 04 August 2005 13:40, you wrote:
It might be good to have an example of your problem. On Aug 4, 2005, at 5:57 AM, Mark Miller wrote:
Hi, I have been trying to figure out how to use the nonlinear regression to fit the cumulative lognormal distribution to a number of data points I have but I am a new R user and I cant quite decipher the notes on nonlinear regression. Any help in this regard will be greatly appreciated, my email address is mmiller at nassp.uct.ac.za
______________________________________________ R-help at stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting- guide.html ------------------------------------------------------------------------ rm(list=ls()) outPut = read.csv("dataOut2.csv") arrive = outPut[1] register = outPut[2] complete = outPut[3] IAT = 0 TTR = 0 TTC = 0 TFR = 0 cnt = 1 for(i in array(2:dim(arrive)[1])) { temp = outPut[i,3]-outPut[i,2] if(temp > 0) { IAT[cnt] = outPut[i,1]-outPut[i-1,1] TTR[cnt] = outPut[i,2]-outPut[i,1] TFR[cnt] = outPut[i,3]-outPut[i,2] cnt = cnt + 1 } } cumIAT = IAT/sum(IAT) for(i in array(2:length(IAT))) { cumIAT[i] = cumIAT[i-1]+cumIAT[i] } cumIAT[1] = 0 plot(cumIAT,do.point=F) TTR[cnt] = outPut[1,2]-outPut[1,1] TFR[cnt] = outPut[1,3]-outPut[1,2] # Plot for inter-arrival times x = seq(0,30,0.01) #postscript("cumIAT.ps") plot(ecdf(IAT), do.point=FALSE) lines(x, pexp(x,0.4)) dev.off() # rexp(100,0.21) x = seq(0,20,0.01) postscript("cumTTR.ps") plot(ecdf(TTR), do.point=FALSE) lines(x, plnorm(x,1,0.7)) dev.off() # rlnorm(100,1,0.7) # Plot for Time to complete from registered x = seq(0,30,0.01) postscript("cumTFR.ps") plot(ecdf(TFR), do.point=FALSE) lines(x*600, pbeta(x,1.4,4.3)) dev.off() # rbeta(100,1.6,5)*600 # Find the position with the leat time and hence the next avaliable ambulance minimum = function(toFind) { min = 0; pos = 0; for(i in array(1:length(toFind))) { if(i == 1) { min = toFind[i] pos = i } else { if(toFind[i] < min) { min = toFind[i] pos = i } } } pos } ambsReq = 0 numAmbs = 0 numberAmbs = 0 avgWait = 1 numberAmbs2 = 0 avgWaitTime2 = 0 avgWaitTime = 0 counter = 0 counter2 = 1 cntO = 1 for(i in array(1:50)) { while(avgWait > 0) { counter = counter + 1 numAmbs = numAmbs + 1 numberAmbs[counter] = numAmbs numCalls = 1 ambs = array(c(array(0,numAmbs)), dim=c(numAmbs,numCalls)) waitTime = ambs totalTime = ambs currTime = 0 timeTS = 0 IotherAT = 0 TotherTR = 0 for(i in array(1:500)) { #interAT = IAT(ceil(rand()*length(IAT))); interAT = rexp(1,0.21) #timeTR = TTR(ceil(rand()*length(TTR))); timeTR = rlnorm(1,1,0.7) #timeFR = TFR(ceil(rand()*length(TFR))); timeFR = rbeta(1,1.4,4.3)*600 IotherAT[i] = interAT TotherTR[i] = timeTR currTime = currTime + interAT pos = minimum(totalTime) if(ambs[pos,numCalls] != 0) { numCalls = numCalls + 1 ambs = array(c(ambs,array(0,numAmbs)), dim=c(numAmbs,numCalls)) waitTime = array(c(waitTime,array(0,numAmbs)), dim=c(numAmbs,numCalls)) if(totalTime[pos] > currTime) waitTime[pos,numCalls] = totalTime[pos] - currTime totalTime[pos] = totalTime[pos] + interAT + timeTR + timeFR ambs[pos,numCalls] = interAT + timeTR + timeFR } else { for(i in array(1:numCalls)) { if(ambs[pos,i] == 0) { if(totalTime[pos] > currTime) waitTime[pos,i] = totalTime[pos] - currTime totalTime[pos] = totalTime[pos] + interAT + timeTR + timeFR ambs[pos,numCalls] = interAT + timeTR + timeFR break } } } } avgWait = sum(waitTime)/500 avgWaitTime[counter] = avgWait if(numAmbs == 25) { avgWaitTime2[cntO] = avgWait numberAmbs2[cntO] = numAmbs cntO = cntO + 1 } } postscript("timeAmbs.ps") plot(numberAmbs,avgWaitTime,'lines') dev.off() postscript("timeAmbs2.ps") plot(numberAmbs2,avgWaitTime2,'lines') dev.off() ambsReq[counter2] = numAmbs numAmbs = 0 numberAmbs = 0 avgWait = 1 avgWaitTime = 0 counter = 0 counter2 = counter2 +1 cntO = 1 numberAmbs2 = 0 avgWaitTime2 = 1 } TotherFR = 0 cnt = 1 for(i in array(1:dim(ambs)[1])) { for(j in array(1:dim(ambs)[2])) { if(ambs[i,j] != 0) { TotherFR[cnt] = ambs[i,j] cnt = cnt +1 } } } postscript("dataIAT.ps") hist(IAT,30) dev.off() postscript("simIAT.ps") hist(IotherAT,30) dev.off() postscript("dataTTR.ps") hist(TTR,30) dev.off() postscript("simTTR.ps") hist(TotherTR,30) dev.off() postscript("dataTFR.ps") hist(TFR,30) dev.off() postscript("simTFR.ps") hist(TotherFR,30) dev.off() ------------------------------------------------------------------------ ______________________________________________ R-help at stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Spencer Graves, PhD Senior Development Engineer PDF Solutions, Inc. 333 West San Carlos Street Suite 700 San Jose, CA 95110, USA spencer.graves at pdf.com www.pdf.com <http://www.pdf.com> Tel: 408-938-4420 Fax: 408-280-7915