Estimate cross-covariance of seasonality (NDVI) and precipitation data with a harmonic regression model then visualize a chart from the model

The challenge is to estimate cross-covariance of seasonality (NDVI) and 30-day precipitation with a harmonic model. I added the UCSB CHIRPS rainfall data as bands to all available Landsat 8 imagery and computed cross-covariance. Would it be possible to plot the covariance with the original NDVI data and obtain a chart of the plot. I have attempted this but received an error. The step at Line 183 creates an image and later when I try to map a function over that image it doesn’t work. I believe .map() will only work with a Feature Collection or Image Collection.

How do I get the covariance as a Image Collection, then use it in the harmonic model?

I have tried this code:

// var collection = ee.ImageCollection.fromImages(image.get('images'));
// return ee.Image(image).addBands(collection.mean().rename('mean'));

I have included the script.

Can you give me advice on what I might have to add or change?

Script is here: https://code.earthengine.google.com/ee49496c6bbc2ab3a83898c348d0ae48