headers
Justus Schwabedal | 1 Apr 02:07
Picon
Picon

Re: Estimation of parameters while fitting data

The Problem you are dealing with is not easily solved, especially,  
when you are dealing with an nonlinear differential equation. However  
there are methods available though it involves a little background  
knowledge. I can give you some literature:

Braake Braake Bock Briggs: Fitting ordinary differential equations to  
chaotic data, Phys. Rev. A (1990)

These guys use a special kind of least squares fitting, namely a  
multiple shooting method. In their work it seems to work reasonable,  
however I think it is very unstable under the influence of any noise.

Parlitz Junge Kocarev: Synchronization-based estimation of parameters  
from time series, PRE (1996)

I wrote my diploma thesis on this method and were able to estimate  
parameters from EEG time series of epileptic patients. The idea is,  
that the model will synchronize with the data when the model  
parameters fit the (suspected) parameters of the system most well. If  
your model fits the data well you should stick to this. I would send  
you my thesis but it's written in german so it probably won't be of  
any help.
If you have questions, please ask. Yours Justus

On Mar 31, 2008, at 11:51 PM, Anne Archibald wrote:

> On 31/03/2008, Doreen Mbabazi <doreen <at> aims.ac.za> wrote:
>
>> Thanks, I tried to do that(by taking err = V-f(y,t,p)[2]) while  
>> defining
(Continue reading)

Permalink | Reply | 0 comments |
headers
Perry Greenfield | 1 Apr 03:59

Re: FITS images with header-supplied axes?


On Mar 30, 2008, at 12:38 PM, Keflavich wrote:
>
> I was thinking no resampling, just put an RA/DEC grid and fit the
> image into it as well as it can be.  I don't know if it's possible to
> display rotated pixels, but that would be the most useful behavior in
> this case.
>
Just to clarify, do you mean nearest neighbor regridding? It does  
sound like you mean that you would like to see the image rotated to  
have north up. Right?

> I'm not sure I understand how the transforms would make it easier to
> display sky coordinates without labeling axes, though.  Are you saying
> that if the image is already sampled in RA/DEC (or whatever coordinate
> system) space, then it should be easy to display the RA/DEC
> coordinates on the axes?
>
The transforms machinery would allow displaying the image in its  
original orientation and then when the cursor was moved over the  
image, the displayed x, y coordinates could display RA, DEC instead  
of pixel coordinates. But if you desire to rotate the image and have  
it aligned with north, that capability isn't really important. Doing  
the rotated image needs some tool to do the rotation (quickly  
presumably rather than a precise mapping) and then display that as  
part of a larger tool. That would be much more straightforward. But  
we haven't made such a tool yet. If someone does it first, it sure  
would be nice to add, at least as part of an astronomy toolkit.

Perry
(Continue reading)

Permalink | Reply | 3 comments |
headers
Robert Kern | 1 Apr 05:17
Picon
Gravatar

Re: kmeans2 random initialization

On Mon, Mar 31, 2008 at 1:19 PM, Alec Solway <asolway <at> seas.upenn.edu> wrote:
> Hello,
>
>  In using kmeans2 with the "random" initialization method, I often get:
>
>  numpy.linalg.linalg.LinAlgError: Matrix is not positive definite -
>      Cholesky decomposition cannot be computed
>
>  I haven't dug too deeply into the implementation; I'm wondering if
>  there's a gotcha someone can point out quickly here. What exactly is
>  the "random" initialization method doing with the data set, and what
>  constraints must the data have so as to keep whatever intermediate
>  matrix it computes to be positive definite?

The relevant function is scipy/cluster/vq.py:_krandinit(). It is
finding the covariance matrix and manually doing a multivariate normal
sampling. Your data is most likely degenerate and not of full rank.
It's arguable whether or not this should fail, but
numpy.random.multivariate_normal() uses the SVD instead of a Cholesky
decomposition to find the matrix square root, so it sort of ignores
non-positive definiteness.

Try replacing the last 2 lines of _krandinit() with

  x = N.random.multivariate_normal(mu, cov, k)

and see if that helps you.

--

-- 
Robert Kern
(Continue reading)

Permalink | Reply | 2 comments |
headers
David Warde-Farley | 1 Apr 10:41
Picon
Favicon
Gravatar

Re: kmeans2 random initialization

On 31-Mar-08, at 11:17 PM, Robert Kern wrote:

> The relevant function is scipy/cluster/vq.py:_krandinit(). It is
> finding the covariance matrix and manually doing a multivariate normal
> sampling. Your data is most likely degenerate and not of full rank.
> It's arguable whether or not this should fail, but
> numpy.random.multivariate_normal() uses the SVD instead of a Cholesky
> decomposition to find the matrix square root, so it sort of ignores
> non-positive definiteness.

This might not be relevant, depending on how the covariance is  
computed, but one 'gotcha' I've seen with numerical algorithms that  
assume positive-definiteness is that occasionally floating point  
oddities will induce (very slight) non-symmetry of the input matrix,  
and thus the algorithm will choke; it's easily solved by averaging the  
matrix with it's transpose (though there are probably more efficient  
ways).

David
Permalink | Reply | 1 comment |
headers
humufr | 1 Apr 16:41
Picon
Favicon

Re: FITS images with header-supplied axes?

				Hello,

I just want to discuss about a problem that our python astronomer have and the 
answer from Perry is a very good example of it. There are too many projects 
which are doing exactly the same thing and essentially because nothing is 
centralized or worst, in the case of pywcs, not advertised. 

For example Perry told us about pywcs developped by STSCI (it was the first 
time I saw any reference to this project) but Adam spoke about 
astlib ( http://astlib.sourceforge.net/ ) which are different package with 
exactly the same goal and I did myself something similar (even if it was fast 
and dirtier).

I think it's time to try to identified the most important task that the 
astronomer are needing and try to centralized all the effort at the same 
place. The astropy mail list is probably a good start as the 
astropy.scipy.org website.

Perhaps we can start to identified the need and desiderata. We will have too  
many but that some most important or urgent must be identified. As example, 
we clearly needed the pywcs and thanks to STSCI we have it now. We also need 
a package to plot our data, images etc. Matplotlib is very good but does have 
a major problem (at least for my point of view), it's slow, very slow for big 
array or can't even produce the image if the image is too big. ( If I 
remember a precedent discussion, the problem is mainly due to Agg.)

After we need to know, if possible, for which project if someone is doing 
something, who and how to contact him/her. So interested people can 
eventually help him/her to do it. It's seems that STSCI is doing most of the 
work but I'm pretty sure that other people are willingly to help them to
(Continue reading)

Permalink | Reply | 0 comments |
headers
heytrent@gmail.com | 1 Apr 18:36
Picon

Real-time plotting and data storage format questions

Greetings,

A small team of us are developing a new simulation package from the ground up. Our legacy approach relied on MATLAB and other proprietary software. A hope pf ours is to be able to shed the use of MATLAB for the analysis of our simulation results and instead use python with scipy/numpy/matplotlib etc. I've successfully installed and compiled optimized numpy/scipy and all the supporting packages (ATLAS, FFTW, etc).

So far so good.

To the point - I have two questions:

1) We would like to have a "scope" to monitor simulation outputs in real time. We're using one tool that can take data over a tcp/ip port, but is clunky and only works on a single platform. Does such a thing exists within the python realm for plotting data in real time?

2) Our simulation creates large (1-4 GB) data sets. Since we're writing this simulation ourselves (C++) we can save the data in any format. Does anyone have a suggestion for a specific format or API that's been found to be optimal in terms of memory usage and ability to import into python for analysis and plotting?

Thank you for any suggestions. We're still new with Python, so I apologize if these questions seem mundane.

_______________________________________________
SciPy-user mailing list
SciPy-user <at> scipy.org
http://projects.scipy.org/mailman/listinfo/scipy-user
Permalink | Reply | 14 comments |
headers
Alan G Isaac | 1 Apr 19:01
Picon
Favicon
Gravatar

Re: Real-time plotting and data storage format questions

1. Will this project produce code that you
will be sharing?  If so, please post updates
to this list!

2.  Is the following of interest?
Last example at 
<URL:http://code.enthought.com/chaco/gallery/index.shtml>
Code at
<URL:http://pyserial.sourceforge.net/>
Comments at:
<URL:https://svn.enthought.com/enthought/browser/Chaco/trunk/examples/advanced/spectrum.py>

Cheers,
Alan Isaac
Permalink | Reply | 0 comments |
headers
Charles Doutriaux | 1 Apr 19:06

Re: Real-time plotting and data storage format questions

Hi there,

For the data format I strongly recommend NetCDF and especially NetCDF4 
(allows for compression)
In addition of NetCDF4 I strongly recommend you look at the CF 
conventions, which is a standard and allows software to learn about the 
data you're storing, see:
http://cf-pcmdi.llnl.gov/

As far as an all in 1 python based software to store/analyze/plot your 
data, you can look at CDAT (which happens to redistribute scipy )

Hope this helps,

C.

heytrent <at> gmail.com wrote:
> Greetings,
>
> A small team of us are developing a new simulation package from the 
> ground up. Our legacy approach relied on MATLAB and other proprietary 
> software. A hope pf ours is to be able to shed the use of MATLAB for 
> the analysis of our simulation results and instead use python with 
> scipy/numpy/matplotlib etc. I've successfully installed and compiled 
> optimized numpy/scipy and all the supporting packages (ATLAS, FFTW, etc).
>
> So far so good.
>
> To the point - I have two questions:
>
> 1) We would like to have a "scope" to monitor simulation outputs in 
> real time. We're using one tool that can take data over a tcp/ip port, 
> but is clunky and only works on a single platform. Does such a thing 
> exists within the python realm for plotting data in real time?
>
> 2) Our simulation creates large (1-4 GB) data sets. Since we're 
> writing this simulation ourselves (C++) we can save the data in any 
> format. Does anyone have a suggestion for a specific format or API 
> that's been found to be optimal in terms of memory usage and ability 
> to import into python for analysis and plotting?
>
> Thank you for any suggestions. We're still new with Python, so I 
> apologize if these questions seem mundane.
> ------------------------------------------------------------------------
>
> _______________________________________________
> SciPy-user mailing list
> SciPy-user <at> scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-user
>
Permalink | Reply | 10 comments |
headers
Gael Varoquaux | 1 Apr 19:09
Favicon
Gravatar

Re: Real-time plotting and data storage format questions

On Tue, Apr 01, 2008 at 10:06:35AM -0700, Charles Doutriaux wrote:
> For the data format I strongly recommend NetCDF and especially NetCDF4 
> (allows for compression)

Actually, why not hdf5, which seems to be used by NetCDF4, but is also 
largely used across many scientific comunities, and very well supported
under python (pytables)?

My 2 cents,

Gaël
Permalink | Reply | 9 comments |
headers
Ryan Krauss | 1 Apr 20:16
Picon

Re: Impulse/Step response troubles

I have seen the problem with repeated roots.  I thought I reported it.

Your second case with a ramp response could be better handled with
signal.lsim using your original transfer function ([K],[1.0,8.0,K])
and u = t.

On Mon, Mar 31, 2008 at 11:50 AM, Kumar Appaiah <akumar <at> iitm.ac.in> wrote:
> Reply below quote.
>
>
>  On Mon, Mar 31, 2008 at 08:52:57AM +0530, Kumar Appaiah wrote:
>  > The question is to determine (plot) the step and ramp responses of
>  > K / (s^2 + 8s + K) for K = 7, 16 and 80
>  >
>  > So, the following code is used, and I'll keep changing a line:
>  > <code>
>  > from scipy import *
>  > from pylab import *
>  >
>  > K = 16.0 # should be redone with 7.0, 16.0 and 80.0
>  >
>  > r = signal.impulse(([K],[1.0,8.0,K]), T=r_[0:5:0.001])
>  > plot(r[0], r[1], linewidth=2)
>  > show()
>  > </code>
>  >
>  > The above code plots the impulse response of the given function. Now,
>  > on to the action:
>  >
>  > 1. Running the above code with K = 7.0 and K = 16.0 gives expected
>  > results. However, running the code with K = 16.0 doesn't work; if K =
>  > 16.0, the response should be 16 * t * exp(-4t) * u(t), which it
>  > isn't. Changing K to 16 + or - .00000000001 fixes it. There surely is
>  > some problem when the value of K is such that a double root is hit.
>
>  Well, I've zeroed in on the issue. Running the code of
>  signal.lti.impulse, we get here:
>
>  s,v = linalg.eig(sys.A)
>  vi = linalg.inv(v)
>
>  Now, v is a matrix which has the eigen vectors, which are EQUAL in
>  this case:
>
>  print sys.A
>  [[-0.9701425  -0.9701425 ]
>   [ 0.24253563  0.24253563]]
>
>  which is singular. However, the code goes on to invert this happily,
>  resulting in a bad matrix and horrendous values. I would be surprised
>  if nobody has encountered this yet (didn't find this on the Tickets).
>
>  What would be the best way to handle repeated roots in the transfer
>  function's denominator?
>
>  Thanks.
>
>
>
>  Kumar
>  --
>  Kumar Appaiah,
>  458, Jamuna Hostel,
>  Indian Institute of Technology Madras,
>  Chennai - 600 036
>  _______________________________________________
>  SciPy-user mailing list
>  SciPy-user <at> scipy.org
>  http://projects.scipy.org/mailman/listinfo/scipy-user
>
Permalink | Reply | 1 comment |

Gmane