[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