2

I have the following formula:

enter image description here

that i have to code (C is a mtrix, y and y_star are vectors). I tried the following, according to this post: Python numpy array sum over certain indices

> def chisq02(y, y_star, C_inv):
>     indices = len(y)
>     return np.sum(np.sum(np.matmul(y[i] - y_star[i], np.matmul(C_inv, y[j]-y_star[j])) for i in indices)for j in indices)

but it returns me the error:

TypeError: 'int' object is not iterable

Obviously it does not work... Does anyone know how this could work ?

Ehsan
  • 12,072
  • 2
  • 20
  • 33
Apinorr
  • 166
  • 3
  • 14

3 Answers3

2

Your indices variable is just an integer, the length of y. But it must be a list from 0 to len(y) to loop over all indices. You can do this like so:

indices_list = range(len(y))

The range() function produces a list [0, 1, ..., len(y)-1]. Then you can loop over this list in the list comprehension.

Bram Dekker
  • 631
  • 1
  • 7
  • 19
1

indices = len(y) returns you an integer, so to iterate you need to use for i in range(indices).

I suspect you don't need to use np.matmul in your case because you do scalar multiplication according to your formula. I suppose C_inv[i][j] should be used instead of C_inv.

Oleh
  • 91
  • 3
  • Yes i think you're right ! I tried this, and it didn't return me an error: np.sum( np.sum(y[i] - y_star[i]* C_inv[i][j]* y[j]-y_star[j] for i in indices)for j in indices) – Apinorr Jul 02 '20 at 10:46
1

If your matrices/vectors are not arrays, first convert them to arrays:

import numpy as np

y = np.array(y)
y_star = np.array(y_star)
C_inv = np.array(C_inv)

Use this, no need for loops:

def chisq02(y, y_star, C_inv):
    return np.sum(np.matmul(np.matmul(y-y_star, C_inv), (y-y_star)))

If your y and y_star are vectors, simply this will work:

def chisq02(y, y_star, C_inv):
    return np.matmul(np.matmul(y-y_star, C_inv), (y-y_star))
Ehsan
  • 12,072
  • 2
  • 20
  • 33
  • But i need to sum two times. Are you sure that this is equivalent to 2 sums ? – Apinorr Jul 02 '20 at 10:41
  • @user11659631 Always try with a sample toy to check. `np.sum` sums all elements by default: checkout numpy doc https://numpy.org/doc/stable/reference/generated/numpy.sum.html: _The default, axis=None, will sum all of the elements of the input array._ – Ehsan Jul 02 '20 at 10:43
  • 1
    @user11659631 in fact, you dont even need sum. just remove it. – Ehsan Jul 02 '20 at 10:45
  • I am so confused, why don't I need to sum ?:/ – Apinorr Jul 02 '20 at 10:51
  • 1
    Bacause the matrix multiplication `np.matmul` has the sum inside it. When you multiply matrices, you multiply rows by columns and then sum. That is the whole point of using numpy and matrix multiplication. looping with numpy is discouraged and generally not a good idea. If you have lists convert them to numpy first. I will add the edit to post for it. – Ehsan Jul 02 '20 at 10:53
  • And matrix multiplication always start from the rigth to the left so should it be np.matmul(y-y_star, np.matmul(C_inv, y-y_star)) – Apinorr Jul 02 '20 at 10:53
  • Ok i see I think i need to think that trough ! Thank you ! – Apinorr Jul 02 '20 at 10:55
  • 1
    matrix multiplication has associative property, it does not matter which you do first: A(BC) = (AB)C – Ehsan Jul 02 '20 at 10:55