Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I just want to perform the following operation done in Matlab into the Field Calculator using Python:

p(p==0) = 0.001 ; p(p == 1) = 0.999; 
tVar = .5*log(p ./ (1 - p))

the operation consists of a Logit transformation:

  1. Find 0's and 1's and replace with .001 and .999 respectively, and
  2. Perform the second calculation with logarithms

I've tried the following for 1) in Python to no avail:

Pre-Logic Script Code:

def iLogit(a):
  return [.0001 if x==0 else x for x in a]

Any help please?

share|improve this question
    
A better formula for this purpose--because it would be a monotonic transformation of the data, which yours is not--would simply add a small constant to each of p and 1-p, as in 0.5 * log((p + 0.001) / (1 - p + 0.001)). If your value of p is found as a ratio of counts, you would do even better by recomputing those ratios after adding a small "start" value (such as 1/6) to both the numerator and the denominator. This will prevent p from ever attaining 0 or 1. –  whuber Jan 12 at 18:00

3 Answers 3

up vote 6 down vote accepted

Using a Python parser:

Pre-Logic Script Code:

def switch(x):
    if x == 0:
        x = 0.001
    elif x == 1:
        x = 0.999
    tVar = 0.5 * math.log(x / (1 - x))
    return tVar

In the Codeblock:

y=

switch(!x!)

enter image description here

share|improve this answer
1  
There is no need for the else block if you do nothing. –  Nathan W Jan 12 at 23:41
    
Thanks @Nathan, I made that change. Simple is better than complex. –  Aaron Jan 13 at 1:51

Try with...

Pre-Logic Script Code:

def iLogit(x):
  if x==0:
    x=0.001
  elif x==1:
    x=0.999
  tVar=0.5*math.log(x/(1-x))
  return tVar

EDIT (compact version):

def iLogit(x):
    x = .001 if x==0 else .999 if x==1 else x
    return .5*math.log(x/(1-x))
share|improve this answer
    
Thanks a lot @afalciano! –  Oliver Amundsen Jan 29 at 11:40

updated answer

As you use a list (return [.0001 if x==0 else x for x in a]) and a a logit funtion (thanks whuber, I had completely forgotten this function in my first answer: real values strictly between 0 and 1.):

in Matlab/Octave:

 p=[1,1,0.5,0.6,0,0.4,1];
 p(p==0) = 0.001 ; p(p == 1) = 0.999;
 disp(p)
 0.9990    0.9990    0.5000    0.6000    0.0010    0.4000    0.9990 
 tVar = .5*log(p ./ (1 - p));
 disp(tvar)
 3.4534    3.4534         0    0.2027   -3.4534   -0.2027    3.4534

In Python:

from math import log #natural logarithm
p=[1,1,0.5,0.6,0,0.4,1]
# with conditional list comprehension
p =[(0.001 if x==0 else (0.999 if x==1  else x)) for x in p]
print p
[0.999, 0.999, 0.5, 0.6, 0.001, 0.4, 0.999]
tVar = [.5*log(x/(1-x)) for x in p]
print tVar
[3.4533773893242765, 3.4533773893242765, 0.0, 0.2027325540540821, -3.453377389324277, -0.20273255405408214, 3.4533773893242765]

So, we can write the function (result as a list):

def iLogit(a):
    replace = [(0.001 if x==0 else (0.9999 if x==1  else x)) for x in a]
    return [.5*log(x/(1-x)) for x in replace]

( Original answer with errors in the values for logit transformation )

in Matlab:

p=[1,2,0,1,0,2,1,1];
p(p==0) = 0.001 ; p(p == 1) = 0.999;
disp(p):
0.9990000   2.0000000   0.0010000   0.9990000   0.0010000   2.0000000 0.9990000   0.9990000
tVar = .5*log(p ./ (1 - p));
disp(tvar)
3.45338 + 0.00000i   0.34657 + 1.57080i  -3.45338 + 0.00000i   3.45338 + 0.00000i  -3.45338 + 0.00000i   0.34657 + 1.57080i   3.45338 + 0.00000i   3.45338 + 0.00000i

In Python:

from math import log #natural logarithm
p=[1,2,0,1,0,2,1,1]
# with conditional list comprehension
p =[(0.001 if x==0 else (0.999 if x==1  else x)) for x in p]
print p
[0.999, 2, 0.001, 0.999, 0.001, 2, 0.999, 0.999]
tVar = [.5*log(x/(1-x)) for x in p]

So, in theory:

def iLogit(a):
    replace = [(0.001 if x==0 else (0.9999 if x==1  else x)) for x in a]
    return [.5*log(x/(1-x)) for x in replace]

Why in theory ? Because of the logarithm of negative numbers: you can only compute the logarithm of a positive number in Python only.

 x = 3
 print .5*log(x/(1-x)) # or print .5*log(3/(1-3))
 Traceback (most recent call last):
     File "<stdin>", line 1, in <module>
 ValueError: math domain error

 # with numpy
 print .5*numpy.log(x/(1-x))
 Warning: invalid value encountered in log
 nan

If you want the same result as Matlab (complex number), you must use the cmath module

import cmath
x = 0.999
print .5*(cmath.log(x/(1-x)))
(3.45337738932+0j)

Don't forget the problem of floating Point Arithmetic: What Every Computer Scientist Should Know About Floating-Point Arithmetic and Floating Point Arithmetic: Issues and Limitations in Python.)

share|improve this answer
    
A logit transformation is intended to be applied only to real values strictly between 0 and 1. –  whuber Jan 12 at 21:18
    
Sorry, I don't know ArcPy, so my answer. Then my function iLogit(a):works without problem. –  gene Jan 12 at 21:26
    
Whether this is implemented in ArcPy or not is irrelevant. See en.wikipedia.org/wiki/Logit, for instance. –  whuber Jan 12 at 21:39
    
Thanks, whuber, I had completely forgotten this function... –  gene Jan 13 at 8:34
    
Amazing explanation Thanks! @gene –  Oliver Amundsen Jan 29 at 11:40

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.