4

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?

1
  • 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 '14 at 18:00
6

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

2
  • 1
    There is no need for the else block if you do nothing. – Nathan W Jan 12 '14 at 23:41
  • Thanks @Nathan, I made that change. Simple is better than complex. – Aaron Jan 13 '14 at 1:51
3

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.)

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

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))
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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