[R] Using nonlinear regression
Spencer Graves
spencer.graves at pdf.com
Mon Aug 15 17:37:52 CEST 2005
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
More information about the R-help
mailing list