› › ›
Options

Blog - exploring climate data (part 1)

John BaezJohn Baez
in Azimuth Blog

I want us to create some short blog articles that publicize the work you guys have been doing. The idea is to post graphics and descriptions of how they were created, not necessarily drawing big conclusions.

To get started, let's try this:

Dara: could you

  1. Give me the figures that are included in this file?

  2. Say where you got the Darwin and Tahiti data?

  3. Say anything you can to explain how these figures were generated? (I understand how the Gabor transform works - is that all you are using here, or something more complicated? A Gabor transform needs a certain parameter that describes the width of the "window" - what are you using there?)

Ideally people will learn to write these posts themselves; right now we're limited by the fact that I'm doing most of the writing. But I'll try to write this one. There's nothing here yet, but this is where I'll write it:

«1

Comments

  • Options
    1.
    Dara O Shayda
    edited July 2014

    Hello John

    Files

    I added the actual tabular data + in the readme file sepcified a bit more info about where it came from and how it was processed.

    GaborWavelet [4] (order 4) no padding GaborWavelet

    DGaussianWavelet[4] (order 4) no padding DGaussianWavelet

    The transform:

    ContinuousWaveletTransform

    Scalogram (2D plots) Scalogram

    The time-series were transformed:

    DaubechiesWavelet[4] order 4, Padding -> "Extrapolated" DaubechiesWavelet

    DiscreteWaveletTransform

    Let me know if you need more, soon I post the code into github.

    Dara

    Comment Source:Hello John [Files](http://files.lossofgenerality.com/Baez.zip) I added the actual tabular data + in the readme file sepcified a bit more info about where it came from and how it was processed. GaborWavelet [4] (order 4) no padding [GaborWavelet](http://reference.wolfram.com/language/ref/GaborWavelet.html) DGaussianWavelet[4] (order 4) no padding [DGaussianWavelet](http://reference.wolfram.com/language/ref/DGaussianWavelet.html) The transform: [ContinuousWaveletTransform](http://reference.wolfram.com/language/ref/ContinuousWaveletTransform.html) Scalogram (2D plots) [Scalogram](http://reference.wolfram.com/language/ref/WaveletScalogram.html) The time-series were transformed: DaubechiesWavelet[4] order 4, Padding -> "Extrapolated" [DaubechiesWavelet](http://reference.wolfram.com/language/ref/DaubechiesWavelet.html) [DiscreteWaveletTransform](http://reference.wolfram.com/language/ref/DiscreteWaveletTransform.html) Let me know if you need more, soon I post the code into github. Dara
  • Options
    2.
    Dara O Shayda
    edited July 2014

    Scalograms are like piano-roll sheet-music, therefore I like to investigate a dynamical system's data not by human intuition and cognitive observations e.g. 12 month periodicity known to all of us, rather by looking at the scalograms which are basically the time-frequency map of a 1D time-series.

    These are difficult to obtain from FFT and therefore as motivation the new field of Wavelet Theory was developed. FFT basis are cos and sin like with no compact support, but the wavelets' compact support allow for such Scalograms. We will use this property later on for forecast algorithms via the Kernel Mathematics (I will publish a tutorial very soon).

    Also I do multiple decompositions (since they are approximate in nature) to render a clearer understanding of the time-frequency map.

    I have converted the frequency to periods (months) i.e. the y-axis is periods in months.

    Dara

    Comment Source:Scalograms are like piano-roll sheet-music, therefore I like to investigate a dynamical system's data not by human intuition and cognitive observations e.g. 12 month periodicity known to all of us, rather by looking at the scalograms which are basically the time-frequency map of a 1D time-series. These are difficult to obtain from FFT and therefore as motivation the new field of Wavelet Theory was developed. FFT basis are cos and sin like with no compact support, but the wavelets' compact support allow for such Scalograms. We will use this property later on for forecast algorithms via the Kernel Mathematics (I will publish a tutorial very soon). Also I do multiple decompositions (since they are approximate in nature) to render a clearer understanding of the time-frequency map. I have converted the frequency to periods (months) i.e. the y-axis is periods in months. Dara
  • Options
    3.
    Dara O Shayda

    John wrote:

    not necessarily drawing big conclusions.

    To be on the clear, I am not interested in the conclusion or inference parts of this endeavor. Strides in computing models and robust well back-tested forecast algorithms are top priority in my mind + I love to learn from John as much as possible.

    Dara

    Comment Source:John wrote: >not necessarily drawing big conclusions. To be on the clear, I am not interested in the conclusion or inference parts of this endeavor. Strides in computing models and robust well back-tested forecast algorithms are top priority in my mind + I love to learn from John as much as possible. Dara
  • Options
    4.
    John BaezJohn Baez
    edited July 2014

    Thanks, Dara! I'll try to write a blog article based on this, and I'll try very hard to keep it short and simple (I'm afraid that may be hard).

    I'll let everyone know here when it's ready for your improvements!

    Comment Source:Thanks, Dara! I'll try to write a blog article based on this, and I'll try very hard to keep it short and simple (I'm afraid that may be hard). I'll let everyone know here when it's ready for your improvements!
  • Options
    5.
    John BaezJohn Baez

    Okay, I'm ready for comments here:

    Some comments of my own:

    1. I do not understand the formula for the wavelet transforms being used here. If I could write down the mathematical formula I would understand what's going on. Mathematica code would be helpful too, though I'd still need to translate it into equations.

    2. The use of a "frequency 4" Gabor wavelet and the 4th derivative of a Gaussian seems arbitrary and mysterious to me.

    3. Why is there almost no power at periods shorter than 9.5 months? Was the data "denoised" before feeding it into these wavelet transforms? If so, how? Does that explain the shortage of power at low periods? I would expect a lot of power at periods all the way down to 1 month.

    4. I understand the windowed Fourier transform or Gabor transform much better than these various wavelet transforms. I would love to see a Gabor transform of this Darwin and Tahiti transform. I believe the Gabor transform is approximately the same as a continuous wavelet transform using a "frequency 1" Gabor wavelet. I'm not sure.

    Comment Source:Okay, I'm ready for comments here: * [[Blog - exploring climate data (part 1)]] Some comments of my own: 1. I do not understand the formula for the wavelet transforms being used here. If I could write down the mathematical formula I would understand what's going on. Mathematica code would be helpful too, though I'd still need to translate it into equations. 1. The use of a "frequency 4" Gabor wavelet and the 4th derivative of a Gaussian seems arbitrary and mysterious to me. 1. Why is there almost no power at periods shorter than 9.5 months? **Was the data "denoised" before feeding it into these wavelet transforms? If so, how?** Does that explain the shortage of power at low periods? I would expect a lot of power at periods all the way down to 1 month. 1. I understand the windowed Fourier transform or [Gabor transform](http://johncarlosbaez.wordpress.com/2013/01/30/milankovich-vs-the-ice-ages/) much better than these various wavelet transforms. I would love to see a Gabor transform of this Darwin and Tahiti transform. I believe the Gabor transform is approximately the same as a continuous wavelet transform using a "frequency 1" Gabor wavelet. I'm not sure.
  • Options
    6.
    Graham Jones

    Just a typo:

    Puzzle 1. Why are there more low-frequency (i.e., long-period) fluctuations in Tahiti and Darwin?

    I guess you mean "Tahiti than Darwin"

    Comment Source:Just a typo: > Puzzle 1. Why are there more low-frequency (i.e., long-period) fluctuations in Tahiti and Darwin? I guess you mean "Tahiti than Darwin"
  • Options
    7.
    DavidTweed
    edited July 2014

    HI John, looking at the wiki page it seems to me that the "inline images" for the 4th derivatives of Gaussian don't match the linked images (which I looked at to get a bigger look). The inline ones look like they might be exactly the same as the Gabor wavelet images, but I'm not confident enough in this kind of wavelet analysis to say that's not right :-)

    EDIT: I actually looked at the html source, saw the paste problem and took the liberty of fixing it.

    A general question: the vertical axis of the graphs is a log-scale. Does that affect how I should interpret the graph? (I don't know, I'm asking.)

    Comment Source:HI John, looking at the wiki page it seems to me that the "inline images" for the 4th derivatives of Gaussian don't match the linked images (which I looked at to get a bigger look). The inline ones look like they might be exactly the same as the Gabor wavelet images, but I'm not confident enough in this kind of wavelet analysis to say that's not right :-) EDIT: I actually looked at the html source, saw the paste problem and took the liberty of fixing it. A general question: the vertical axis of the graphs is a log-scale. Does that affect how I should interpret the graph? (I don't know, I'm asking.)
  • Options
    8.
    WebHubTelWebHubTel

    At some point we should note that the Tahiti-Darwin pair was selected because it approximates a dipole, which is the embodiment of a perfect -1 anti-correlation that you mention.

    Yet, Darwin and Tahiti were also chosen because they have a long history of almost continuous atmospheric pressure readings, and actually don't represent the best dipole pair in the Pacific. I have seen Port Moresby mentioned as a better candidate than Darwin, while Tahiti is not as good as it can be either.

    I am working extra-curricularly with the climatologists and computer scientists at the U of Minnesota (a very interesting mix) and they have some interesting ideas on how to find the ideal dipoles from the geospatial and temporal records:

    [1]J. Kawale, S. Liess, A. Kumar, M. Steinbach, A. R. Ganguly, N. F. Samatova, F. H. Semazzi, P. K. Snyder, and V. Kumar, “Data Guided Discovery of Dynamic Climate Dipoles.,” presented at the CIDU, 2011, pp. 30–44.

    http://www.researchgate.net/publication/228459734_Data_Guided_Discovery_of_Dynamic_Climate_Dipoles/file/72e7e5182f2be780a0.pdf

    Comment Source:At some point we should note that the Tahiti-Darwin pair was selected because it approximates a dipole, which is the embodiment of a perfect -1 anti-correlation that you mention. Yet, Darwin and Tahiti were also chosen because they have a long history of almost continuous atmospheric pressure readings, and actually don't represent the best dipole pair in the Pacific. I have seen Port Moresby mentioned as a better candidate than Darwin, while Tahiti is not as good as it can be either. I am working extra-curricularly with the climatologists and computer scientists at the U of Minnesota (a very interesting mix) and they have some interesting ideas on how to find the ideal dipoles from the geospatial and temporal records: [1]J. Kawale, S. Liess, A. Kumar, M. Steinbach, A. R. Ganguly, N. F. Samatova, F. H. Semazzi, P. K. Snyder, and V. Kumar, “Data Guided Discovery of Dynamic Climate Dipoles.,” presented at the CIDU, 2011, pp. 30–44. <http://www.researchgate.net/publication/228459734_Data_Guided_Discovery_of_Dynamic_Climate_Dipoles/file/72e7e5182f2be780a0.pdf>
  • Options
    9.
    Nathan Urban

    Cool, I attended Kumar's talk on their dipole work at the 2012 Climate Informatics workshop. Are you extending this work specifically?

    Comment Source:Cool, I attended Kumar's talk on their dipole work at the 2012 Climate Informatics workshop. Are you extending this work specifically?
  • Options
    10.
    Dara O Shayda
    edited July 2014

    Hello John

    With regards to your questions, first draft of some detailed answers:

    Sanity Check for Wavelets

    i. The last plots in the pdf are the denoised data which was not included in the original stuff I sent you. ii. if you need the images from this file, all are available in jpg could zip and send you

    General thoughts:

    1. I like to hammer the code and keep looking for bugs, the more we question the code the better, happier I am :)
    2. I love writing live code tutorials to make sure I learn myself and could explain the concepts, so there is always ancillary coding for that matter. There are geniuses out there who just know things, I am not one of them, so when I talk about code it is about learning and aiding in thinking about the research, not about a final thing with definite veracity i.e. lots of bugs ;)
    3. Wavelets like anything else have artifacts and they need to be learned and experimented with e.g. aliasing and padding.
    4. Change of definitions: People keep talking about periodicity of El Nino for example, I understand Local Periodicity in order to be able to deal with discrete finite historical data. Inferences must be made cognizant of the fact that frequency or periodicity or amplitude are localized in time (or the varying parameter of the t-series function). At such time t, forecast-able periodicity is x , if you remove the localization then I do not understand the arguments and inferences

    I go back to your blog and try to work on your questions and concerns.

    Dara

    Comment Source:Hello John With regards to your questions, first draft of some detailed answers: [Sanity Check for Wavelets](http://files.lossofgenerality.com/sanityWLT.pdf) i. The last plots in the pdf are the denoised data which was not included in the original stuff I sent you. ii. if you need the images from this file, all are available in jpg could zip and send you General thoughts: 1. I like to hammer the code and keep looking for bugs, the more we question the code the better, happier I am :) 2. I love writing live code tutorials to make sure I learn myself and could explain the concepts, so there is always ancillary coding for that matter. There are geniuses out there who just know things, I am not one of them, so when I talk about code it is about learning and aiding in thinking about the research, not about a final thing with definite veracity i.e. lots of bugs ;) 3. Wavelets like anything else have artifacts and they need to be learned and experimented with e.g. aliasing and padding. 4. Change of definitions: People keep talking about periodicity of El Nino for example, I understand Local Periodicity in order to be able to deal with discrete finite historical data. Inferences must be made cognizant of the fact that frequency or periodicity or amplitude are localized in time (or the varying parameter of the t-series function). At such time t, forecast-able periodicity is x , if you remove the localization then I do not understand the arguments and inferences I go back to your blog and try to work on your questions and concerns. Dara
  • Options
    11.
    John BaezJohn Baez
    edited July 2014

    Another question, Dara. Exactly what time period do you use for your Darwin and Tahiti data? I'm a bit worried because of how the red line shoots down at the right end here:

    If you look at the air pressure data files you sent me, the Darwin file ends like this:

    2012   6.0   7.2   6.3  10.7  11.3  13.3  12.8  14.0  12.2  10.8   9.3   7.7
2013   5.8   7.2   7.2   9.7  10.8  11.2  12.5  13.1  11.6  10.7   7.7   7.7
2014   5.0 -99.9 -99.9 -99.9 -99.9 -99.9 -99.9 -99.9 -99.9 -99.9 -99.9 -99.9
2015 -99.9 -99.9 -99.9 -99.9 -99.9 -99.9 -99.9 -99.9 -99.9 -99.9 -99.9 -99.9

    while Tahiti ends like this:

    2012  12.4  12.3  10.8  12.2  12.6  13.2  13.5  14.9  15.0  14.1  12.8  10.0
2013  10.0  11.0  13.2  12.1  13.5  14.1  14.7  14.7  14.6  13.3 -99.9 -99.9

    The number -99.9 means "no data". But if you don't stop your data analysis in October 2013, you'll get some large negative values for the air pressure!

    Comment Source:Another question, Dara. Exactly what time period do you use for your Darwin and Tahiti data? I'm a bit worried because of how the red line shoots down at the right end here: <img src = "http://math.ucr.edu/home/baez/ecological/shayda/darwin_tahiti_old.jpg" alt = ""/> If you look at the air pressure data files you sent me, the Darwin file ends like this: ~~~~ 2012 6.0 7.2 6.3 10.7 11.3 13.3 12.8 14.0 12.2 10.8 9.3 7.7 2013 5.8 7.2 7.2 9.7 10.8 11.2 12.5 13.1 11.6 10.7 7.7 7.7 2014 5.0 -99.9 -99.9 -99.9 -99.9 -99.9 -99.9 -99.9 -99.9 -99.9 -99.9 -99.9 2015 -99.9 -99.9 -99.9 -99.9 -99.9 -99.9 -99.9 -99.9 -99.9 -99.9 -99.9 -99.9 ~~~~ while Tahiti ends like this: ~~~~ 2012 12.4 12.3 10.8 12.2 12.6 13.2 13.5 14.9 15.0 14.1 12.8 10.0 2013 10.0 11.0 13.2 12.1 13.5 14.1 14.7 14.7 14.6 13.3 -99.9 -99.9 ~~~~ The number -99.9 means "no data". But if you don't stop your data analysis in October 2013, you'll get some large negative values for the air pressure!
  • Options
    12.
    John BaezJohn Baez
    edited July 2014

    Dara wrote:

    first draft of some detailed answers:

    Sanity Check for Wavelets

    Unfortunately this pdf doesn't answer my questions. Here are my questions, phrased more clearly I hope. Without answers to these it'll be hard for me to finish this blog article. Can you please answer them here on the forum?

    1. What is the formula for the continuous wavelet transforms being used here?

    2. What is a "frequency 4" Gabor wavelet?

    3. Why are you using a frequency 4 Gabor wavelet?

    4. Why are you using the 4th derivative of a Gaussian?

    5. Why is there almost no power at periods shorter than 9.5 months?

    6. Did you "denoise" the time series before feeding them into these wavelet transforms?

    7. If so, how?

    I think I can answer question 1 using the Mathematica page Continuous Wavelet Transform. Namely:

    The continuous wavelet transform of a sequence of numbers $x_1, x_2, \dots , x_n$ is

    $$ w(u,s) = \frac{1}{\sqrt{s}} \sum_{k = 1}^n x_k \psi(\frac{\Delta (k-u)}{s}) $$ Here we think of $x_k$ as the value of some function $x(t)$ defined at all times, via

    $$ x_k = x(t_0 + k \Delta) $$ The function $\psi$ is something we get to choose. We could choose it to be a "Gabor wavelet of frequency 4" - but I don't know what this function is, because it's not explained on the Mathematica page. Or we could choose it to the 4th derivative of a Gaussian, as explained here.

    (In all these links you need to click "Details" to get the information I need.)

    The idea behind the continuous wavelet transform is evident to me from the formula above, so if that's what you're using I can explain it. However, I don't see why you chose the frequency 4 Gabor wavelet and the 4th derivative of a Gaussian. These are functions with about 4 wiggles in them; I see why you want some wiggles, but not why you want 4. If there's nothing special about "4", fine, just say so.

    Comment Source:Dara wrote: > first draft of some detailed answers: > [Sanity Check for Wavelets](http://files.lossofgenerality.com/sanityWLT.pdf) Unfortunately this pdf doesn't answer my questions. Here are my questions, phrased more clearly I hope. Without answers to these it'll be hard for me to finish this blog article. Can you please answer them here on the forum? 1. What is the formula for the continuous wavelet transforms being used here? 1. What is a "frequency 4" Gabor wavelet? 1. Why are you using a frequency 4 Gabor wavelet? 1. Why are you using the 4th derivative of a Gaussian? 1. Why is there almost no power at periods shorter than 9.5 months? 1. Did you "denoise" the time series before feeding them into these wavelet transforms? 1. If so, how? I think I can answer question 1 using the Mathematica page [Continuous Wavelet Transform](http://reference.wolfram.com/language/ref/ContinuousWaveletTransform.html). Namely: The continuous wavelet transform of a sequence of numbers $x_1, x_2, \dots , x_n$ is $$ w(u,s) = \frac{1}{\sqrt{s}} \sum_{k = 1}^n x_k \psi(\frac{\Delta (k-u)}{s}) $$ Here we think of $x_k$ as the value of some function $x(t)$ defined at all times, via $$ x_k = x(t_0 + k \Delta) $$ The function $\psi$ is something we get to choose. We could choose it to be a "Gabor wavelet of frequency 4" - but I don't know what this function is, because it's not explained [on the Mathematica page](http://reference.wolfram.com/language/ref/GaborWavelet.html). Or we could choose it to the 4th derivative of a Gaussian, as explained [here](http://reference.wolfram.com/language/ref/GaborWavelet.html). (In all these links you need to click "Details" to get the information I need.) The _idea_ behind the continuous wavelet transform is evident to me from the formula above, so if that's what you're using I can explain it. However, I don't see why you chose the frequency 4 Gabor wavelet and the 4th derivative of a Gaussian. These are functions with about 4 wiggles in them; I see why you want some wiggles, but not why you want 4. If there's nothing special about "4", fine, just say so.
  • Options
    13.
    John BaezJohn Baez

    Graham wrote:

    I guess you mean “Tahiti than Darwin".

    Thanks, I'll fix that.

    Comment Source:Graham wrote: > I guess you mean “Tahiti than Darwin". Thanks, I'll fix that.
  • Options
    14.
    John BaezJohn Baez
    edited July 2014

    I wrote:

    I believe the Gabor transform is approximately the same as a continuous wavelet transform using a “frequency 1” Gabor wavelet. I’m not sure.

    I no longer believe this! If you compare my formula for the continuous wavelet transform to the formula for the Gabor transform you'll see they're very different. And I'm not talking about the fact that the first has a sum while the second has an integral - that's not so important!

    I'm talking about the fact that in a continuous wavelet transform we're convolving our function by a rescaled version of a chosen "wavelet", while in the Gabor transform we're multiplying our function by a fixed Gaussian and then taking its Fourier transform. Even if we choose our wavelet to be a Gaussian, these are different things to do.

    The continuous wavelet transform sees how much a small segment of our function "matches" a rescaled version of our chosen wavelet $\psi$. So, if $\psi$ has some wiggles, this transform will see how much a small segment of our function has wiggles of a given shape and frequency.

    Comment Source:I wrote: > I believe the Gabor transform is approximately the same as a continuous wavelet transform using a “frequency 1” Gabor wavelet. I’m not sure. I no longer believe this! If you compare my formula for the continuous wavelet transform to the formula for the [Gabor transform](https://en.wikipedia.org/wiki/Gabor_transform) you'll see they're very different. And I'm not talking about the fact that the first has a sum while the second has an integral - that's not so important! I'm talking about the fact that in a continuous wavelet transform we're convolving our function by a rescaled version of a chosen "wavelet", while in the Gabor transform we're multiplying our function by a fixed Gaussian and then taking its Fourier transform. Even if we choose our wavelet to be a Gaussian, these are different things to do. The continuous wavelet transform sees how much a small segment of our function "matches" a rescaled version of our chosen wavelet $\psi$. So, if $\psi$ has some wiggles, this transform will see how much a small segment of our function has wiggles of a given shape and frequency.
  • Options
    15.
    Dara O Shayda

    Hello John, no I do not use the undefined, I dropped the 4 rows of data for darwin.ascii and last 2 rows for tahiti.ascii. So you are correct that dangling red tale was caused by the drop last 4 for darwin, I corrected it, only dropping the last 2 rows and putting it here:

    Correction: Tahiti Darwin plots

    Please note that the plots are translated along the y-axis and the plot has aspect ratio of 1/5 or the data is impossible to plot.

    Both data were denoised DaubechiesWavelet[4], 4 levels of refinement, Padding -> "Extrapolated"

    Dara

    Comment Source:Hello John, no I do not use the undefined, I dropped the 4 rows of data for darwin.ascii and last 2 rows for tahiti.ascii. So you are correct that dangling red tale was caused by the drop last 4 for darwin, I corrected it, only dropping the last 2 rows and putting it here: [Correction: Tahiti Darwin plots](http://files.lossofgenerality.com/darwin_tahiti.jpg) Please note that the plots are translated along the y-axis and the plot has aspect ratio of 1/5 or the data is impossible to plot. Both data were denoised DaubechiesWavelet[4], 4 levels of refinement, Padding -> "Extrapolated" Dara
  • Options
    16.
    Dara O Shayda
    edited July 2014

    Please kindly refresh the link in #16 for Correction to get the latest plot.

    I had to cut all the way to 2012, so the dimensions of the data i.e. the number of rows match for Correlation:

    Darwin {2012, 6., 7.2, 6.3, 10.7, 11.3, 13.3, 12.8, 14., 12.2, 10.8, 9.3, 7.7}

    Tahiti {2012, 12.4, 12.3, 10.8, 12.2, 12.6, 13.2, 13.5, 14.9, 15., 14.1, 12.8, 10.}

    The correlation between the raw data then is 0.583649

    Comment Source:Please kindly refresh the link in #16 for Correction to get the latest plot. I had to cut all the way to 2012, so the dimensions of the data i.e. the number of rows match for Correlation: Darwin {2012, 6., 7.2, 6.3, 10.7, 11.3, 13.3, 12.8, 14., 12.2, 10.8, 9.3, 7.7} Tahiti {2012, 12.4, 12.3, 10.8, 12.2, 12.6, 13.2, 13.5, 14.9, 15., 14.1, 12.8, 10.} The correlation between the raw data then is 0.583649
  • Options
    17.
    Dara O Shayda
    edited July 2014

    John wrote:

    The continuous wavelet transform sees how much a small segment of our function “matches” a rescaled version of our chosen wavelet ψ. So, if ψ has some wiggles, this transform will see how much a small segment of our function has wiggles of a given shape and frequency.

    This is what I refer to as the compact support which could be re-scaled.

    Comment Source:John wrote: >The continuous wavelet transform sees how much a small segment of our function “matches” a rescaled version of our chosen wavelet ψ. So, if ψ has some wiggles, this transform will see how much a small segment of our function has wiggles of a given shape and frequency. This is what I refer to as the compact support which could be re-scaled.
  • Options
    18.
    John BaezJohn Baez
    edited July 2014

    Dara wrote:

    So you are correct that dangling red tale was caused by the drop last 4 for darwin, I corrected it, only dropping the last 2 rows...

    Thanks! The new improved graph is here:

    I'm proud that I could guess the old one was a bit wrong.

    Comment Source:Dara wrote: > So you are correct that dangling red tale was caused by the drop last 4 for darwin, I corrected it, only dropping the last 2 rows... Thanks! The new improved graph is here: <img src = "http://math.ucr.edu/home/baez/ecological/shayda/darwin_tahiti.jpg" alt = ""/> I'm proud that I could guess the old one was a bit wrong. <img src = "http://math.ucr.edu/home/baez/emoticons/tongue2.gif" alt = ""/>
  • Options
    19.
    John BaezJohn Baez

    Dara wrote:

    This is what I refer to as the compact support which could be re-scaled.

    In mathematics, a function $f : \mathbb{R} \to \mathbb{R}$ has compact support if it's zero outside some finite interval $[-M,M]$. The Gabor and differentiated Gaussian wavelets we've been discussing now don't have compact support, though they do decrease very rapidly.

    Comment Source:Dara wrote: > This is what I refer to as the compact support which could be re-scaled. In mathematics, a function $f : \mathbb{R} \to \mathbb{R}$ has **compact support** if it's zero outside some finite interval $[-M,M]$. The Gabor and differentiated Gaussian wavelets we've been discussing now don't have compact support, though they do decrease very rapidly.
  • Options
    20.
    John BaezJohn Baez

    Dara wrote:

    The correlation between the raw data then is 0.583649

    Do you mean -0.583649? The normalized cross-correlation should be negative here.

    Comment Source:Dara wrote: > The correlation between the raw data then is 0.583649 Do you mean -0.583649? The <a href = "https://en.wikipedia.org/wiki/Cross-correlation#Normalized_cross-correlation">normalized cross-correlation</a> should be negative here.
  • Options
    21.
    DavidTweed

    I'll just mention what I wasn't sure about on the plots, and might be something other people ask. So the vertical axis of "Gabor/D4OG feature" length is on a log scale. So if I see a pixel of the same colour at one height in the graph and the same colour lower down on the graph, what does that mean? Is it that the pixel describes a sample point at some precise value for the feature length, or is it (an approximation to) the total of the responses for all the feature period values that fall within that pixel. Additionally as the feature length increases its energy will generally increase just because there's more points, unless it's normalised in some way?

    Basically, whether a given colour means the same power wherever it is in the plot?

    Comment Source:I'll just mention what I wasn't sure about on the plots, and might be something other people ask. So the vertical axis of "Gabor/D4OG feature" length is on a log scale. So if I see a pixel of the same colour at one height in the graph and the same colour lower down on the graph, what does that mean? Is it that the pixel describes a sample point at some precise value for the feature length, or is it (an approximation to) the total of the responses for all the feature period values that fall within that pixel. Additionally as the feature length increases its energy will generally increase just because there's more points, unless it's normalised in some way? Basically, whether a given colour means the same power wherever it is in the plot?
  • Options
    22.
    Dara O Shayda

    John wrote:

    Do you mean -0.583649?

    According to Mathematica's Correlation function:

    Correlation

    1. Raw data correlation between Tahiti and Darwin data is 0.583649
    2. But if both denoised by Daubechies [4] index {0,0,0,0} Correlation is -0.507138 a negative number

    3. But if both denoised by Daubechies [4] index {0,0,1} Correlation is 0.911759 i.e. 12 month periodic part of the data

      I checked the code and seems to be in order!

    Dara

    Comment Source:John wrote: >Do you mean -0.583649? According to Mathematica's Correlation function: [Correlation](http://reference.wolfram.com/language/ref/Correlation.html) 1. Raw data correlation between Tahiti and Darwin data is 0.583649 2. But if both denoised by Daubechies [4] index {0,0,0,0} Correlation is -0.507138 a negative number 3. But if both denoised by Daubechies [4] index {0,0,1} Correlation is 0.911759 i.e. 12 month periodic part of the data I checked the code and seems to be in order! Dara
  • Options
    23.
    Dara O Shayda

    David Tweed:

    Basically, whether a given colour means the same power wherever it is in the plot?

    I do not use Log, Mathematica's scalograms are by default in Log/Octave scales. I dislike that, so I wrote a converter to get rid of the log so I could report periodicity with the exact values, here is the code I use:

    period = 1/(cwd[ "SampleRate"]/(#* cwd["Wavelet"]["FourierFactor"])) & /@ (Thread[{Range[ cwd["Octaves"]], 1}] /. cwd["Scales"]); ticks = Transpose[{Range[Length[period]], period}];

    Comment Source:David Tweed: >Basically, whether a given colour means the same power wherever it is in the plot? I do not use Log, Mathematica's scalograms are by default in Log/Octave scales. I dislike that, so I wrote a converter to get rid of the log so I could report periodicity with the exact values, here is the code I use: period = 1/(cwd[ "SampleRate"]/(#* cwd["Wavelet"]["FourierFactor"])) & /@ (Thread[{Range[ cwd["Octaves"]], 1}] /. cwd["Scales"]); ticks = Transpose[{Range[Length[period]], period}];
  • Options
    24.
    Dara O Shayda
    edited July 2014

    John wrote:

    The new improved graph is here:

    John thank you for paying attention to these things and I need that sorta focus on details. However that being said, the plots are padded:

    Array Padding

    This is so since the mathematics needs powers of 2 size data (plus 1?) and most data are not.

    So I use the extrapolated for 1-D time series (I could change that),for 2D scalograms I use NONE, you have these options for {a,b,c} as sample data:

    "Periodic": {c, a, b, c, a, b, c, a, b, c, a}

    "Reversed": {c, c, b, a, a, b, c, c, b, a, a}

    "ReversedNegation": {c, -c, -b, -a, a, b, c, -c, -b, -a, a}

    "Reflected": {a, b, c, b, a, b, c, b, a, b, c}

    "ReflectedDifferences": {2 a + b - 2 c, 2 a - c, 2 a - b, a, b, c, -b + 2 c, -a + 2 c, -2 a + b + 2 c}

    "Extrapolated" {4 a - 3 b, 3 a - 2 b, 2 a - b, a, b, c, -b + 2 c, -2 b + 3 c, -3 b + 4 c}

    Comment Source:John wrote: > The new improved graph is here: John thank you for paying attention to these things and I need that sorta focus on details. However that being said, the plots are padded: [Array Padding](http://reference.wolfram.com/language/ref/ArrayPad.html) This is so since the mathematics needs powers of 2 size data (plus 1?) and most data are not. So I use the extrapolated for 1-D time series (I could change that),for 2D scalograms I use NONE, you have these options for {a,b,c} as sample data: "Periodic": {c, a, b, c, a, b, c, a, b, c, a} "Reversed": {c, c, b, a, a, b, c, c, b, a, a} "ReversedNegation": {c, -c, -b, -a, a, b, c, -c, -b, -a, a} "Reflected": {a, b, c, b, a, b, c, b, a, b, c} "ReflectedDifferences": {2 a + b - 2 c, 2 a - c, 2 a - b, a, b, c, -b + 2 c, -a + 2 c, -2 a + b + 2 c} "Extrapolated" {4 a - 3 b, 3 a - 2 b, 2 a - b, a, b, c, -b + 2 c, -2 b + 3 c, -3 b + 4 c}
  • Options
    25.
    Dara O Shayda

    John wrote:

    The Gabor and differentiated Gaussian wavelets we’ve been discussing now don’t have compact support, though they do decrease very rapidly.

    Sigh... thanx again, people use this terms for compact support as t goes to infinity and often loosely... I do not know if there is a better way of saying this.

    What I mean is a burst short-lived periodic signal

    Dara

    Comment Source:John wrote: >The Gabor and differentiated Gaussian wavelets we’ve been discussing now don’t have compact support, though they do decrease very rapidly. Sigh... thanx again, people use this terms for compact support as t goes to infinity and often loosely... I do not know if there is a better way of saying this. What I mean is a burst short-lived periodic signal Dara
  • Options
    26.
    Dara O Shayda

    John

    Correction to the time-series data did not change the scalograms noticeably but I will reissue and post them again tonight.

    Dara

    Comment Source:John Correction to the time-series data did not change the scalograms noticeably but I will reissue and post them again tonight. Dara
  • Options
    27.
    Dara O Shayda

    Hello John

    if you need more compuations or plots let me know, this stuff was just to start something, but let me know what you think is necessary for an investigation

    Dara

    Comment Source:Hello John if you need more compuations or plots let me know, this stuff was just to start something, but let me know what you think is necessary for an investigation Dara
  • Options
    28.
    John BaezJohn Baez
    edited July 2014

    Dara wrote:

    According to Mathematica’s Correlation function:

    Correlation

    Okay, the "details" here say Correlation computes Pearson's correlation coefficient, which seem to be exactly the same as what I was calling the normalized cross-correlation.

    Namely, for two lists of numbers $v_i$, $w_i$ with the same length, this is

    $$ \frac{\langle v w \rangle \;-\; \langle v \rangle \langle w \rangle}{sd(v) \; sd(w) } $$ where $\langle \rangle $ means 'mean' and $sd$ means 'standard deviation'.

    Raw data correlation between Tahiti and Darwin data is 0.583649

    That's impossible with the definition I just gave! You can just look: the correlation has to be negative, because the blue curve here tends to be above zero when the red one is below zero, and vice versa:

    So, there is some mistake somewhere.

    Comment Source:Dara wrote: > According to Mathematica’s Correlation function: > [Correlation](http://reference.wolfram.com/language/ref/Correlation.html) Okay, the "details" here say Correlation computes [Pearson's correlation coefficient](https://en.wikipedia.org/wiki/Pearson%27s_correlation_coefficient#Definition), which seem to be exactly the same as what I was calling the [normalized cross-correlation](https://en.wikipedia.org/wiki/Cross-correlation#Normalized_cross-correlation). Namely, for two lists of numbers $v_i$, $w_i$ with the same length, this is $$ \frac{\langle v w \rangle \;-\; \langle v \rangle \langle w \rangle}{sd(v) \; sd(w) } $$ where $\langle \rangle $ means 'mean' and $sd$ means 'standard deviation'. > Raw data correlation between Tahiti and Darwin data is 0.583649 That's impossible with the definition I just gave! You can just look: the correlation has to be negative, because the blue curve here tends to be above zero when the red one is below zero, and vice versa: <img src = "http://math.ucr.edu/home/baez/ecological/shayda/darwin_tahiti.jpg" alt = ""/> So, there is some mistake somewhere.
  • Options
    29.
    John BaezJohn Baez

    Dara, I would be much happier if you'd just go through my list of questions and answer them all. "I don't know" is an acceptable answer. However, without knowing the answers to these, I don't know what you're doing or why you're doing it:

    1. What is the formula for the continuous wavelet transforms being used here?

    2. What is a "frequency 4" Gabor wavelet?

    3. Why are you using a frequency 4 Gabor wavelet?

    4. Why are you using the 4th derivative of a Gaussian?

    5. Why is there almost no power at periods shorter than 9.5 months?

    6. Did you "denoise" the time series before feeding them into these wavelet transforms?

    7. If so, how?

    I guessed an answer to 1; you haven't commented on that guess.

    You've answered 6 and 7.

    Comment Source:Dara, I would be much happier if you'd just go through my list of questions and answer them all. "I don't know" is an acceptable answer. However, without knowing the answers to these, I don't know what you're doing or why you're doing it: 1. What is the formula for the continuous wavelet transforms being used here? 1. What is a "frequency 4" Gabor wavelet? 1. Why are you using a frequency 4 Gabor wavelet? 1. Why are you using the 4th derivative of a Gaussian? 1. Why is there almost no power at periods shorter than 9.5 months? 1. Did you "denoise" the time series before feeding them into these wavelet transforms? 1. If so, how? I guessed an answer to 1; you haven't commented on that guess. You've answered 6 and 7.
  • Options
    30.
    Dara O Shayda
    edited July 2014

    John could you go to these links:

    1.What is the formula for the continuous wavelet transforms being used here? ContinuousWaveletTransform

    And then click on DETAILS AND OPTIONS look for square bullet item "The continuous wavelet transform of a function..." and the formula is there. Sorry cannot typeset here easily.

    2.What is a “frequency 4” Gabor wavelet? Gabor Wavelet

    click on DETAILS second square bullet is the formula, replace w by 4, if you need that done symbolically in a few moments I post a little thing after this...

    3.Why are you using a frequency 4 Gabor wavelet? I chose 4 for no particular reason other than it gives a non-even Scalogram, if you increase the frequency to 6 and more then solid frequency bands appear and I feel that it might be more revealing if a lesser frequency used. This part of choosing the parameters for wavelets is mostly experimentation, I did not find any references on how to figure these numbers out algebraically. I do not make any assumptions just try different numbers and report the results or match to what is known already.

    As a matter of fact I plot and analyse using many frequencies, but I did not want to clutter the documents, by all means one should use many frequencies and compare results. If I do it my way, I will generate a large set of such plots in the server then someone like yourself cruise through and decide which frequency is best to consider.

    4.Why are you using the 4th derivative of a Gaussian? Same as 3, but usually I like to use 40 and more! to get a better sense of the burst-noise in the data.

    5.Why is there almost no power at periods shorter than 9.5 months?

    Example power number: 0.0018 or 0.002 for {1} is still significant! usually the NOISE for me is 10^-6 and I consider 10^-3 significant. Therefore somehow I am not sure how to think about this. If you are asking me why these numbers are not larger, I could not give you a single response, I could only say that is what is in the data. Best possible guess the data collection issue.

    Dara

    Comment Source:John could you go to these links: 1.What is the formula for the continuous wavelet transforms being used here? [ContinuousWaveletTransform](http://reference.wolfram.com/language/ref/ContinuousWaveletTransform.html) And then click on DETAILS AND OPTIONS look for square bullet item "The continuous wavelet transform of a function..." and the formula is there. Sorry cannot typeset here easily. 2.What is a “frequency 4” Gabor wavelet? [Gabor Wavelet](http://reference.wolfram.com/language/ref/GaborWavelet.html) click on DETAILS second square bullet is the formula, replace w by 4, if you need that done symbolically in a few moments I post a little thing after this... 3.Why are you using a frequency 4 Gabor wavelet? I chose 4 for no particular reason other than it gives a non-even Scalogram, if you increase the frequency to 6 and more then solid frequency bands appear and I feel that it might be more revealing if a lesser frequency used. This part of choosing the parameters for wavelets is mostly experimentation, I did not find any references on how to figure these numbers out algebraically. I do not make any assumptions just try different numbers and report the results or match to what is known already. As a matter of fact I plot and analyse using many frequencies, but I did not want to clutter the documents, by all means one should use many frequencies and compare results. If I do it my way, I will generate a large set of such plots in the server then someone like yourself cruise through and decide which frequency is best to consider. 4.Why are you using the 4th derivative of a Gaussian? Same as 3, but usually I like to use 40 and more! to get a better sense of the burst-noise in the data. 5.Why is there almost no power at periods shorter than 9.5 months? Example power number: 0.0018 or 0.002 for {1} is still significant! usually the NOISE for me is 10^-6 and I consider 10^-3 significant. Therefore somehow I am not sure how to think about this. If you are asking me why these numbers are not larger, I could not give you a single response, I could only say that is what is in the data. Best possible guess the data collection issue. Dara
  • Options
    31.
    Dara O Shayda

    John

    I plotted the two wavelet families with different frequencies, if you need to use the images I could post the jpg files:

    Gabor and DGaussian frequencies

    Comment Source:John I plotted the two wavelet families with different frequencies, if you need to use the images I could post the jpg files: [Gabor and DGaussian frequencies](http://files.lossofgenerality.com/wavelets.pdf)
  • Options
    32.
    John BaezJohn Baez
    edited July 2014

    Thanks for answering my questions, Dara! I'm a "word guy" - when people ask me questions I like trying to answer them with exactly the right words, and when I ask questions I really want people to do the same. You weren't doing that. Now you did! I think now I can finish the blog article in an hour or two.

    Comment Source:Thanks for answering my questions, Dara! I'm a "word guy" - when people ask me questions I like trying to answer them with exactly the right words, and when I ask questions I really want people to do the same. You weren't doing that. Now you did! I think now I can finish the blog article in an hour or two.
  • Options
    33.
    Dara O Shayda

    You weren’t doing that. Now you did!

    I did not see your questions earlier yesterday, sorry when I code which I do now full-time, I get very scattered and sloppy with anything else that is not code. So do not think I was ignoring you or had an agenda, just did not see or process them in timely manner.

    I’m a “word guy” - when people ask me questions I like trying to answer them with exactly the right words, and when I ask questions I really want people to do the same.

    I must admit this is my issue, I am ok when I am not coding but at this time, please ask me again, normally I am here all the time, if I missed something, I could respond in less than a day.

    Dara

    Comment Source:> You weren’t doing that. Now you did! I did not see your questions earlier yesterday, sorry when I code which I do now full-time, I get very scattered and sloppy with anything else that is not code. So do not think I was ignoring you or had an agenda, just did not see or process them in timely manner. >I’m a “word guy” - when people ask me questions I like trying to answer them with exactly the right words, and when I ask questions I really want people to do the same. I must admit this is my issue, I am ok when I am not coding but at this time, please ask me again, normally I am here all the time, if I missed something, I could respond in less than a day. Dara
  • Options
    34.
    John BaezJohn Baez

    I understand the way it works with coding - not from myself but with other people. My wife used to do software documentation, and it could be hard to get the programmers to explain what they were doing with words. One said "I don't do natural language processing".

    I'm making good progress on the blog article now. Two requests:

    • Could you give me an image of the graph produced by

    Plot[WaveletPsi[DGaussianWavelet[4], x], {x, -5,0,5}]

    I believe this is one of the wavelets you are using.

    • Could you also do the same for the "frequency 4 Gabor wavelet" you are using? Or whatever it is you're using here? I'm a bit confused because this one has both a real and imaginary part. Maybe you are just using the real part?
    Comment Source:I understand the way it works with coding - not from myself but with other people. My wife used to do software documentation, and it could be hard to get the programmers to explain what they were doing with words. One said **"I don't do natural language processing"**. I'm making good progress on the blog article now. Two requests: * Could you give me an image of the graph produced by `Plot[WaveletPsi[DGaussianWavelet[4], x], {x, -5,0,5}]` I believe this is one of the wavelets you are using. * Could you also do the same for the "frequency 4 Gabor wavelet" you are using? Or whatever it is you're using here? I'm a bit confused because this one has both a real and imaginary part. Maybe you are just using the real part?
  • Options
    35.
    Dara O Shayda

    I’m a bit confused because this one has both a real and imaginary part. Maybe you are just using the real part?

    Complex Wavelet calculation

    Let me know if this is enough or I could code something better

    Comment Source:>I’m a bit confused because this one has both a real and imaginary part. Maybe you are just using the real part? [Complex Wavelet calculation ](http://files.lossofgenerality.com/complexWLT.pdf) Let me know if this is enough or I could code something better
  • Options
    36.
    Dara O Shayda

    DGaussian 4

    Blue for Re and yellow/gold for Im Gabor 4

    Comment Source:[DGaussian 4](http://files.lossofgenerality.com/DGauss4.jpg) Blue for Re and yellow/gold for Im [Gabor 4](http://files.lossofgenerality.com/Gabor4.jpg)
  • Options
    37.
    WebHubTelWebHubTel

    Nathan Urban asked:

    "Cool, I attended Kumar’s talk on their dipole work at the 2012 Climate Informatics workshop. Are you extending this work specifically?"

    I was initially trying to get some expert opinion on some questions I had and was directed to the work on dipoles. We are working together in a limited sense and am curious to see where I can take the discussion.

    Comment Source:Nathan Urban asked: "Cool, I attended Kumar’s talk on their dipole work at the 2012 Climate Informatics workshop. Are you extending this work specifically?" I was initially trying to get some expert opinion on some questions I had and was directed to the work on dipoles. We are working together in a limited sense and am curious to see where I can take the discussion.
  • Options
    38.
    John BaezJohn Baez

    Dara wrote:

    Let me know if this is enough or I could code something better.

    Thanks, that's enough! Just right.

    Comment Source:Dara wrote: > Let me know if this is enough or I could code something better. Thanks, that's enough! Just right.
  • Options
    39.
    WebHubTelWebHubTel

    A couple of weeks ago I did the Tahiti-Darwin correlation coefficient calculation that John is talking about here: http://contextearth.com/2014/07/17/correlation-of-time-series/

    The number I get using Mathematica is -0.55, in line with John's interpretation of it being negative negative number because the dipole is anti-correlated.

    In the figure below, I take the negative of the Tahiti measure to make it positive (I also multiply the CC by 100 for visual checking of numbers).

    SOI

    It is possible that Dara is doing something similar by way of taking the complement and just losing track of the sign.

    The number of -0.55 does not look particularly strong and is a weaker anti-correlation than the -0.58, but this number is also extremely sensitive to the amount of high-frequency noise in the time series data. The higher frequencies can rapidly reduce the anti-correlation because any slight amount of phase shifting cause the two waves to rapidly lose alignment according to the CC formula.

    I was hoping to find some other measure that considers the amplitude error caused by slight phase shifts and takes this into account for a "goodness of fit" criteria. In other words the error would be the distance in amplitude and time, instead of just amplitude differences at a particular time. I am not sure that such a beast exists, perhaps related to a variation of a Mahalanobis distance measure? http://en.wikipedia.org/wiki/Mahalanobis_distance

    Comment Source:A couple of weeks ago I did the Tahiti-Darwin correlation coefficient calculation that John is talking about here: <http://contextearth.com/2014/07/17/correlation-of-time-series/> The number I get using Mathematica is -0.55, in line with John's interpretation of it being negative negative number because the dipole is anti-correlated. In the figure below, I take the negative of the Tahiti measure to make it positive (I also multiply the CC by 100 for visual checking of numbers). ![SOI](http://imagizer.imageshack.us/a/img539/1412/48fe30.gif) It is possible that Dara is doing something similar by way of taking the complement and just losing track of the sign. The number of -0.55 does not look particularly strong and is a weaker anti-correlation than the -0.58, but this number is also extremely sensitive to the amount of high-frequency noise in the time series data. The higher frequencies can rapidly reduce the anti-correlation because any slight amount of phase shifting cause the two waves to rapidly lose alignment according to the CC formula. I was hoping to find some other measure that considers the amplitude error caused by slight phase shifts and takes this into account for a "goodness of fit" criteria. In other words the error would be the distance in amplitude and time, instead of just amplitude differences at a particular time. I am not sure that such a beast exists, perhaps related to a variation of a Mahalanobis distance measure? <http://en.wikipedia.org/wiki/Mahalanobis_distance>
  • Options
    40.
    DavidTweed
    edited July 2014

    Hi John, I see that in the blog article you're building up the correlation coefficient from covariance. That's probably the right thing to do, particularly given the other parts of the El Nino series, but just checking that you're aware that you can also develop some ideas about correlation using the following:

    If we've got some time series of data $(x_i)$ and a fixed total amount of "energy" in $y$-stuff equal to the energy in the $x_i$s (say due to normalisation) that we can distribute into a time series $(y_i)$, then it's clear that looking at $\sum_{i=1:N} (x_i-y_i)^2$ gives small values for a very good match and progressively bigger values as they are become less correlated (intuitive meaning of the word). Observing $$ \sum_{i=1:N} (x_i-y_i)^2 =\sum_{i=1:N} (x_i^2 - 2 x_i y_i + y_i^2)=\sum_{i=1:N} x_i^2 - 2 \sum_{i=1:N} x_i y_i + \sum_{i=1:N} y_i^2 $$ Since $\sum_{i=1:N} x_i^2$ is a constant and we've assumed that the total energy $\sum_{i=1:N} y_i^2$ is fixed, we can see that $\sum_{i=1:N} x_i y_i$ contains all the varying behaviour of the $\sum_{i=1:N} (x_i-y_i)^2$ measure, only moving in the opposite direction (ie, big values indicate good matches).

    That's probably actually a more complicated setup than the covariance/correlation one, but I thought I'd mention it.

    Comment Source:Hi John, I see that in the blog article you're building up the correlation coefficient from covariance. That's probably the right thing to do, particularly given the other parts of the El Nino series, but just checking that you're aware that you can also develop some ideas about correlation using the following: If we've got some time series of data $(x_i)$ and a fixed total amount of "energy" in $y$-stuff equal to the energy in the $x_i$s (say due to normalisation) that we can distribute into a time series $(y_i)$, then it's clear that looking at $\sum_{i=1:N} (x_i-y_i)^2$ gives small values for a very good match and progressively bigger values as they are become less correlated (intuitive meaning of the word). Observing \[ \sum_{i=1:N} (x_i-y_i)^2 =\sum_{i=1:N} (x_i^2 - 2 x_i y_i + y_i^2)=\sum_{i=1:N} x_i^2 - 2 \sum_{i=1:N} x_i y_i + \sum_{i=1:N} y_i^2 \] Since $\sum_{i=1:N} x_i^2$ is a constant and we've assumed that the total energy $\sum_{i=1:N} y_i^2$ is fixed, we can see that $\sum_{i=1:N} x_i y_i$ contains all the varying behaviour of the $\sum_{i=1:N} (x_i-y_i)^2$ measure, only moving in the opposite direction (ie, big values indicate good matches). That's probably actually a more complicated setup than the covariance/correlation one, but I thought I'd mention it.
  • Options
    41.
    Graham Jones

    Southern Oscillation Index based upon annual standardization, Climate Analysis Section, NCAR/UCAR. This includes links to monthly sea level pressure anomalies in Darwin and Tahiti, in either ASCII format (click the second two links) or netCDF format (click the first one and read the explanation).

    Which data sets are being used here? There are several at http://www.cgd.ucar.edu/cas/catalog/climind/soiAnnual.html#download. The links in the blog post are to darwin.ascii and tahiti.ascii. From just looking at the numbers, these do not seem to be seasonally adjusted and they do look positively correlated. They both go up around July and down around January.

    Comment Source:> <a href = "http://www.cgd.ucar.edu/cas/catalog/climind/soiAnnual.html">Southern Oscillation Index based upon annual standardization</a>, Climate Analysis Section, NCAR/UCAR. This includes links to monthly sea level pressure anomalies in <a href = "http://www.cgd.ucar.edu/cas/catalog/climind/darwin.ascii">Darwin</a> and <a href = "http://www.cgd.ucar.edu/cas/catalog/climind/tahiti.ascii">Tahiti</a>, in either ASCII format (click the second two links) or netCDF format (click the first one and read the explanation). Which data sets are being used here? There are several at [http://www.cgd.ucar.edu/cas/catalog/climind/soiAnnual.html#download](http://www.cgd.ucar.edu/cas/catalog/climind/soiAnnual.html#download). The links in the blog post are to `darwin.ascii` and `tahiti.ascii`. From just looking at the numbers, these do not seem to be seasonally adjusted and they do look positively correlated. They both go up around July and down around January.
  • Options
    42.
    John BaezJohn Baez

    Graham wrote:

    Which data sets are being used here? There are several at http://www.cgd.ucar.edu/cas/catalog/climind/soiAnnual.html#download. The links in the blog post are to darwin.ascii and tahiti.ascii. From just looking at the numbers, these do not seem to be seasonally adjusted and they do look positively correlated. They both go up around July and down around January.

    Whoops! That would explain it. Those are the data sets Dara is using!

    Sigh. This was supposed to be a quick and easy post, but I guess there's no such thing.

    This would also be consistent with Dara getting huge amounts of "power" in his scalogram at a "period" on the order of one year: he's seeing seasons. ("Period" in quotes because I now realize it's only that up to a constant factor of order roughly 1.)

    Comment Source:Graham wrote: > Which data sets are being used here? There are several at http://www.cgd.ucar.edu/cas/catalog/climind/soiAnnual.html#download. The links in the blog post are to darwin.ascii and tahiti.ascii. From just looking at the numbers, these do not seem to be seasonally adjusted and they do look positively correlated. They both go up around July and down around January. Whoops! That would explain it. Those are the data sets Dara is using! Sigh. This was supposed to be a quick and easy post, but I guess there's no such thing. This would also be consistent with Dara getting huge amounts of "power" in his scalogram at a "period" on the order of one year: he's seeing seasons. ("Period" in quotes because I now realize it's only that up to a constant factor of order roughly 1.)
  • Options
    43.
    John BaezJohn Baez
    edited July 2014

    So, Dara: I'm afraid I'm going to ask you to your scalograms and graph using air pressure anomaly data instead of air pressures. With air pressures, it does not relate to the Southern Oscillation Index in a nice way.

    I'm talking about these:

    "darwin-tahiti.jpg"

    "darwin_dgauss.jpg"

    "tahiti_dgauss.jpg"

    It's fine if you don't do the Gabor wavelets: the Gaussian ones are enough for explaining the idea.

    And I'd be happiest if you did not "denoise" the data. For explanatory purposes it's simpler to work with raw data.

    Darwin air pressure anomalies are here, and Tahiti air pressure anomalies are here.

    Don't forget to cut out the data after 2012!

    Comment Source:So, Dara: I'm afraid I'm going to ask you to your scalograms and graph using air pressure _anomaly_ data instead of air pressures. With air pressures, it does not relate to the Southern Oscillation Index in a nice way. I'm talking about these: "darwin-tahiti.jpg" <img width = "300" src = "http://math.ucr.edu/home/baez/ecological/shayda/darwin_tahiti.jpg" alt = ""/> "darwin_dgauss.jpg" <img width = "300" src = "http://math.ucr.edu/home/baez/ecological/shayda/darwin_dgauss.jpg" alt = ""/> "tahiti_dgauss.jpg" <img width = "300" src = "http://math.ucr.edu/home/baez/ecological/shayda/tahiti_dgauss.jpg" alt = ""/> It's fine if you don't do the Gabor wavelets: the Gaussian ones are enough for explaining the idea. And I'd be happiest if you _did not_ "denoise" the data. For explanatory purposes it's simpler to work with raw data. Darwin air pressure anomalies are [here](http://www.cgd.ucar.edu/cas/catalog/climind/darwin.anom.ascii), and Tahiti air pressure anomalies are [here](http://www.cgd.ucar.edu/cas/catalog/climind/tahiti.anom.ascii). Don't forget to cut out the data after 2012!
  • Options
    44.
    John BaezJohn Baez

    David wrote:

    That’s probably the right thing to do, particularly given the other parts of the El Nino series, but just checking that you’re aware that you can also develop some ideas about correlation using the following...

    Thanks! I'd never really thought of correlations that way, though I know and love the formula you used: it's a Jedi math trick for getting ahold of arbitrary inner products when all you know is inner products of stuff with itself.

    I'll add this approach as an "appendix".

    Comment Source:David wrote: > That’s probably the right thing to do, particularly given the other parts of the El Nino series, but just checking that you’re aware that you can also develop some ideas about correlation using the following... Thanks! I'd never really thought of correlations that way, though I know and love the formula you used: it's a Jedi math trick for getting ahold of arbitrary inner products when all you know is inner products of stuff with _itself_. I'll add this approach as an "appendix".
  • Options
    45.
    Dara O Shayda

    Don’t forget to cut out the data after 2012!

    Hello John I will redo the code with the requirements you specified, sometime today/tonight.

    D

    Comment Source:>Don’t forget to cut out the data after 2012! Hello John I will redo the code with the requirements you specified, sometime today/tonight. D
  • Options
    46.
    DavidTweed
    edited July 2014

    The argument I gave above is a bit awkward just because I've primarily seen it used in finding a noisy version of one signal in another signal, in which case you can argue that the sum of squares for both are constant and of the same magnitude. I didn't want to actually edit the blog post while John's working on it, but I've tried to clear up the logic a bit below that can hopefully just be pasted in if its ok.

    Another way to understand correlation

    David Tweed mentioned another approach from signal processing to understanding the quantity

    $$ \langle x y \rangle = \frac{1}{n} \sum_{i = 1}^n x_i y_i $$ If we've got two lists of data $x$ and $y$ that we want to compare to see if they behave similarly, the first thing we ought to do is multiplicatively scale each one so they're of comparable magnitude. There are various possibilities for assigning a scale, but a reasonable one is to ensure they have equal 'energy'

    $$ \sum_{i=1}^n x_i^2 = \sum_{i=1}^n y_i^2 $$ (This can be achieved by dividing each list by its standard deviation, which is equivalent to what was done in the main derivation above.) Once we've done that then it's clear that looking at

    $$\sum_{i=1}^n (x_i-y_i)^2$$ gives small values when they have a very good match and progressively bigger values as they become less similar. Observe that

    $$ \sum_{i=1}^n (x_i-y_i)^2 =\sum_{i=1}^n (x_i^2 - 2 x_i y_i + y_i^2) $$ $$ =\sum_{i=1}^n x_i^2 - 2 \sum_{i=1}^n x_i y_i + \sum_{i=1}^n y_i^2 $$ Since we've scaled things so that $\sum_{i=1}^n x_i^2$ and $\sum_{i=1}^n y_i^2$ are constants, we can see that when $\sum_{i=1}^n x_i y_i$ becomes bigger,

    $$\sum_{i=1}^n (x_i-y_i)^2$$ becomes smaller. So,

    $$\sum_{i=1}^n x_i y_i$$ serves as a measure of how close the lists are, under these assumptions.

    Comment Source:The argument I gave above is a bit awkward just because I've primarily seen it used in finding a noisy version of one signal in another signal, in which case you can argue that the sum of squares for both are constant and of the same magnitude. I didn't want to actually edit the blog post while John's working on it, but I've tried to clear up the logic a bit below that can hopefully just be pasted in if its ok. ---------- <h3> Another way to understand correlation </h3> David Tweed mentioned another approach from signal processing to understanding the quantity $$ \langle x y \rangle = \frac{1}{n} \sum_{i = 1}^n x_i y_i $$ If we've got two lists of data $x$ and $y$ that we want to compare to see if they behave similarly, the first thing we ought to do is multiplicatively scale each one so they're of comparable magnitude. There are various possibilities for assigning a scale, but a reasonable one is to ensure they have equal 'energy' $$ \sum_{i=1}^n x_i^2 = \sum_{i=1}^n y_i^2 $$ (This can be achieved by dividing each list by its standard deviation, which is equivalent to what was done in the main derivation above.) Once we've done that then it's clear that looking at $$\sum_{i=1}^n (x_i-y_i)^2$$ gives small values when they have a very good match and progressively bigger values as they become less similar. Observe that $$ \sum_{i=1}^n (x_i-y_i)^2 =\sum_{i=1}^n (x_i^2 - 2 x_i y_i + y_i^2) $$ $$ =\sum_{i=1}^n x_i^2 - 2 \sum_{i=1}^n x_i y_i + \sum_{i=1}^n y_i^2 $$ Since we've scaled things so that $\sum_{i=1}^n x_i^2$ and $\sum_{i=1}^n y_i^2$ are constants, we can see that when $\sum_{i=1}^n x_i y_i$ becomes bigger, $$\sum_{i=1}^n (x_i-y_i)^2$$ becomes smaller. So, $$\sum_{i=1}^n x_i y_i$$ serves as a measure of how close the lists are, under these assumptions.
  • Options
    47.
    John BaezJohn Baez
    edited July 2014

    Dara - thanks for offering to redo those 3 images! I'm sorry I didn't catch the problem earlier.

    David - thanks for rewriting that. I cut-and-paste it in. I knew that I'd left it a bit rough regarding the explanation of why we can act like $\sum_i x_i^2 = \sum_i y_^2$ are fixed. This is better.

    WebHubTel - can you remind me of your real name? I want to add information about the correlation of the Darwin and Tahiti air pressure anomalies, and if I use your calculation I'd like to cite you by your real name (if that's okay). Where did you get your data and exactly what time period did you use?

    Everyone - the article is approximately finished; see if you like it:

    I should probably add some references to some 'material for further study' on continuous wavelet transforms. I don't know what would be good.

    Comment Source:Dara - thanks for offering to redo those 3 images! I'm sorry I didn't catch the problem earlier. David - thanks for rewriting that. I cut-and-paste it in. I knew that I'd left it a bit rough regarding the explanation of why we can act like $\sum_i x_i^2 = \sum_i y_^2$ are fixed. This is better. WebHubTel - can you remind me of your real name? <img src = "http://math.ucr.edu/home/baez/emoticons/redface.gif" alt = ""/> I want to add information about the correlation of the Darwin and Tahiti air pressure anomalies, and if I use your calculation I'd like to cite you by your real name (if that's okay). Where did you get your data and exactly what time period did you use? **Everyone** - the article is approximately finished; see if you like it: * [[Blog - exploring climate data (part 1)]] I should probably add some references to some 'material for further study' on continuous wavelet transforms. I don't know what would be good.
  • Options
    48.
    DavidTanzer
    edited July 2014

    Very nice article.

    I changed "La Niño" to "La Niña," in the introduction.

    If we call this wavelet y(t), all that matters is that it's a bump with wiggles on it, and that its mean value is 0

    Perhaps say a bit about why in general the wavelet should have a mean of zero. I see from the wikipedia page continuous wavelet transform that there are conditions on the wavelet to make the original signal recoverable from the transform, and these imply that the mean is zero. But if we're just searching for patterns in the data, it still looks useful to allow "wavelets" such as a positive pulse.

    Comment Source:Very nice article. I changed "La Niño" to "La Niña," in the introduction. > If we call this wavelet y(t), all that matters is that it's a bump with wiggles on it, and that its mean value is 0 Perhaps say a bit about why in general the wavelet should have a mean of zero. I see from the wikipedia page [continuous wavelet transform](http://en.wikipedia.org/wiki/Continuous_wavelet_transform) that there are conditions on the wavelet to make the original signal recoverable from the transform, and these imply that the mean is zero. But if we're just searching for patterns in the data, it still looks useful to allow "wavelets" such as a positive pulse.
  • Options
    49.
    DavidTanzer

    Also consider introducing the term scalogram, and expound it in terms of the pictures you've included.

    Comment Source:Also consider introducing the term scalogram, and expound it in terms of the pictures you've included.
  • Options
    50.
    DavidTanzer

    Perhaps add a couple of sentences explaining the motivation that leads you to the wavelet transform, before proceeding to describe its implementation.

    Something like...we'd like to get a picture of the localized frequency content of the signal at each point in time (and why?), which is something that is addressed by "wavelet analysis."

    Comment Source:Perhaps add a couple of sentences explaining the motivation that leads you to the wavelet transform, before proceeding to describe its implementation. Something like...we'd like to get a picture of the localized frequency content of the signal at each point in time (and why?), which is something that is addressed by "wavelet analysis."
Sign In or Register to comment.