[R] R and MatLab implementations of the same model differs
peter dalgaard
pdalgd at gmail.com
Thu Jul 4 23:25:04 CEST 2013
On Jul 4, 2013, at 20:14 , Berend Hasselman wrote:
>
> On 04-07-2013, at 19:56, peter dalgaard <pdalgd at gmail.com> wrote:
>
>>
>> On Jul 4, 2013, at 19:11 , Berend Hasselman wrote:
>>
>>>
>>> On 04-07-2013, at 18:42, Jannetta Steyn <jannetta at henning.org> wrote:
>>>
>>>> Hi Ben and others
>>>>
>>>> I don't quite know how to explain the "doesn't work" in more detail without
>>>> any visual aid.
>>>
>>> You said that R got into an indefinite loop, whatever that maybe.
>>>
>>>> When you run the two scripts it is easy to see the
>>>> difference. MatLab produces a line on x= -55. This is what I expect - a
>>>> more or less straight line. R on the other hand the result drops from -55
>>>> (the initial value) to -80 and then goes up to -71.37092.
>>>>
>>>> The two scripts have exactly the same equations.
>>>
>>>
>>> I don't think so.
>>> In the R script you have
>>>
>>> init = c(v_axon_AB=-55,mNa_axon_AB=1,hNa_axon_AB=0,mK_axon_AB=1)
>>>
>>> That is not the same as in your Matlab script. To make them the same you should replace the line with
>>>
>>> init = c(v_axon_AB=-55,mNa_axon_AB=0,hNa_axon_AB=1,mK_axon_AB=0)
>>>
>>> Using this line gives quite different results. But not the same as Matlab.
>>>
>>> It seems that you ought to carefully check all your equations.
>>> I don't know enough about Matlab/Octave syntax to determine if the results of the two systems are identical.
>>
>> Also, the line
>>
>> iLeak_axon_AB <- gLeak_axon_AB*(v_axon_AB-ELeak_axon_AB)
>>
>> differs from the Matlab code. Changing it gave me all-NA results, though. And the with() construct assumes that init and parms are named vectors; parms is not, so the values are taken from the global environment. It is not clear whether that makes a difference.
>>
>
> Same NA results here. Making parms a named vector makes no difference but is is certainly necessary to do that.
>
> In that case there is probably an inconsistency in the file simulate.m near the end:
> The last non empty line before the line out = … reads
>
> iLeak_axon = gLeak_axon_AB.*(v_axon_AB-ELeak_axon_AB);
>
> and this doesn't agree with the line
>
> iLeak_axon = ELeak_axon_AB*(v_axon_AB-ELeak_axon_AB);
>
> in the function xdot.
>
> Final question for the OP: if the model is supposed to be dynamic, isn't it suspicious that Matlab gives a constant result?
Also, in this block
gNa_axon_AB=300e-3
gK_axon_AB=52.5-3
gLeak_axon_AB=0.0018e-3
the middle line looks like a typo for 52.5e-3 and the 3rd line looks very low -- googling suggests values like (120, 36, 0.3) mS/cm^3, which would be more consistent with a value around 1e-3.
>
> Berend
>
>>
>>>
>>> Berend
>>>
>>>> I have even named the
>>>> variables the same in the two scripts and copied the equations across to
>>>> make sure I haven't made any typos.
>>>>
>>>> Can one attach images to posts? I'll try. The flatline image is the plot
>>>> from MatLab and the other is the plot from R.
>>>>
>>>> Thanks
>>>> Jannetta
>>>>
>>>>
>>>> On 4 July 2013 16:52, Ben Bolker <bbolker at gmail.com> wrote:
>>>>
>>>>> Berend Hasselman <bhh <at> xs4all.nl> writes:
>>>>>
>>>>>>
>>>>>>
>>>>>> On 04-07-2013, at 17:15, Jannetta Steyn <jannetta <at> henning.org>
>>>>> wrote:
>>>>>>
>>>>>>> Hi folks
>>>>>>>
>>>>>>> I have implemented a model of a neuron using Hodgkin Huxley equations
>>>>> in
>>>>>>> both R and MatLab. My preference is to work with R but R is not giving
>>>>> me
>>>>>>> the correct results. I also can't use ode45 as it just seems to go
>>>>> into an
>>>>>>> indefinite loop. However, the MatLab implementation work fine with
>>>>>>> the same
>>>>>>> equations, parameters and initial values and ode45. Below is my R and
>>>>>>> MatLab implementations.
>>>>>>>
>>>>>>
>>>>>> No problem in running your R file. Have plot.
>>>>>> (Mac mini Core2Duo 2.66Ghz running R 3.0.1 Patched
>>>>>> (2013-06-19 r62992) on Mac OS X 10.8.4: 16.5 seconds.)
>>>>>>
>>>>>> Trying to run your Matlab file in Octave.
>>>>>> No success yet because of unavailable ode45.
>>>>>
>>>>> I'm impressed that you (BH) went to the trouble of checking
>>>>> on this vague "doesn't work" question. If you want to go farther
>>>>> I think you can get ode45 for octave by installing the odepkg
>>>>> package:
>>>>>
>>>>> http://octave.sourceforge.net/odepkg/index.html
>>>>>
>>>>> ______________________________________________
>>>>> R-help at r-project.org mailing list
>>>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>>>> PLEASE do read the posting guide
>>>>> http://www.R-project.org/posting-guide.html
>>>>> and provide commented, minimal, self-contained, reproducible code.
>>>>>
>>>>
>>>>
>>>>
>>>> --
>>>>
>>>> ===================================
>>>> Web site: http://www.jannetta.com
>>>> Email: jannetta at henning.org
>>>> ===================================
>>>> <Rplot01.png><MatPlot01.png>______________________________________________
>>>> R-help at r-project.org mailing list
>>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>>>> and provide commented, minimal, self-contained, reproducible code.
>>>
>>> ______________________________________________
>>> R-help at r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>>> and provide commented, minimal, self-contained, reproducible code.
>>
>> --
>> Peter Dalgaard, Professor,
>> Center for Statistics, Copenhagen Business School
>> Solbjerg Plads 3, 2000 Frederiksberg, Denmark
>> Phone: (+45)38153501
>> Email: pd.mes at cbs.dk Priv: PDalgd at gmail.com
>>
>>
>>
>>
>>
>>
>>
>>
>
--
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd.mes at cbs.dk Priv: PDalgd at gmail.com
More information about the R-help
mailing list