2
\$\begingroup\$

The k-index measures the condition of the magnetosphere. It is usually averaged over three hour, so each day has 8 measurements.

  • The planetary k-index (kp-index) is an average of the measures taken in 15 locations in the world,
  • while the estimated kp-index is an average taken over only 8 locations. The data is available here:

https://www.swpc.noaa.gov/products/planetary-k-index

If you take a look at their data for the last 30 days:

https://services.swpc.noaa.gov/text/daily-geomagnetic-indices.txt

You can see that they give the k-index measured in two different locations (Fredericksburg and College), as well as the estimated kp-index. (You can ignore the columns marked "A".)

A single training point is the 8-value vector of the k-index in Frederickburg on a certain day, and its label is the 8-value vector of the estimated kp-index on that same day. Having trained the algorithm on a large number of sample points, a new query is the k-index in Frederickburg on a different day (with no label), and the algorithm must produce an 8-value vector which is its guess for the estimated kp-index on that day. For example, if the query vector was

2 2 1 2 1 2 1 2

(the k-index vector for Fredericksburg yesterday), then the answer

2.00 0.67 0.67 0.67 0.67 1.33 1.00 1.33

would have zero error, as this is exactly the estimated kp-index for that day. Of course, we can't really expect the algorithm to produce a solution with zero error.

My Intention:

An interesting question is can the k-index in either of the two individual locations be used to predict the estimated kp-index?

My Thoughts:

My thought is that this can be done by training a learning model (for example SVR) on some of the data (where a point is an 8-vector from either Fredericksburg or College representing some day, and its label is an 8-vector that is the estimated kp-index for the say day), and then testing its error on the rest of the data. I am programming in python, I might want to use the sklearn.multioutput.MultiOutputRegressor function.

My answer:

Given dataset.txt

#Importing necessary libraries
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split, GridSearchCV
from sklearn.svm import SVR
from sklearn.multioutput import MultiOutputRegressor
from sklearn.metrics import mean_squared_error, r2_score, mean_absolute_error
from sklearn.preprocessing import StandardScaler
import warnings

#Suppressing warnings for a cleaner 
warnings.filterwarnings("ignore")

# Load the dataset
dataset_path = "dataset.txt"
dataset = pd.read_csv(dataset_path, sep="\t", header=None)

# Split the space-separated data into columns
data = dataset[0].str.split(expand=True)

# Exclude unwanted columns: 4,5,6,7 and the last 5 columns
data = data.drop(columns=[3, 4, 5, 6] + list(range(data.shape[1]-5, data.shape[1])))


# Extracting date, kp-index, and Fredericksburg columns
date_columns = data.iloc[:, :3]
kp_index_columns = data.iloc[:, 3:11]
federickburg_columns = data.iloc[:, 11:19]

# Convert individual date columns to a single date column
date_columns['date'] = date_columns[0] + '-' + date_columns[1] + '-' + date_columns[2]

# Transforming Fredericksburg columns to float type
federickburg_transformed = federickburg_columns.astype(float)

# Defining the features (X) and target (y)
X = federickburg_transformed
y = kp_index_columns.astype(float)

# Splitting the data into training and testing sets (80% training, 20% testing)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)


# Initialize the SVR model with a linear kernel
svr = SVR(kernel='linear', C=1)

# Using MultiOutputRegressor to handle multiple target columns
model = MultiOutputRegressor(svr).fit(X_train, y_train)


# Predict on the test data
y_pred = model.predict(X_test)


# Calculate the Mean Squared Error (MSE) and R^2 score for evaluation
mse = mean_squared_error(y_test, y_pred)
r2 = r2_score(y_test, y_pred)

# Displaying the MSE and R^2 scores 
mse, r2 #(1.6405386761259917, 0.04202137156755004)

# Selecting a random sample from our features
random_index = np.random.randint(0, len(X))
sample_fredericksburg_random = X.iloc[random_index].values.reshape(1, -1)

# Predicting the kp-index for the randomly chosen Federicksburg k-index vector
predicted_kp_index_random = model.predict(sample_fredericksburg_random)

# Retrieving the actual kp-index for the chosen sample
actual_kp_index_random = y.iloc[random_index].values.reshape(1, -1)


# Scaling the features and target data for better performance during grid search
scaler_X = StandardScaler().fit(X_train)
scaler_y = StandardScaler().fit(y_train)
X_train_scaled = scaler_X.transform(X_train)
y_train_scaled = scaler_y.transform(y_train)
X_test_scaled = scaler_X.transform(X_test)
y_test_scaled = scaler_y.transform(y_test)

# Defining the hyperparameter grid for SVR
param_grid = {
 'estimator__C': [0.1, 1, 10, 100],
 'estimator__epsilon': [0.01, 0.1, 1],
 'estimator__kernel': ['linear', 'rbf']
}

# Performing grid search to find the best hyperparameters
grid_search = GridSearchCV(MultiOutputRegressor(SVR()), param_grid, cv=5, scoring='neg_mean_squared_error')
grid_search.fit(X_train_scaled, y_train_scaled)


# Using the best estimator found during grid search
best_svr = grid_search.best_estimator_

# Training the model with the best hyperparameters
best_svr.fit(X_train_scaled, y_train_scaled)


# Making predictions on the scaled test data
y_pred_scaled = best_svr.predict(X_test_scaled)

# Transforming the predicted data back to the original scale
y_pred_best = scaler_y.inverse_transform(y_pred_scaled)

# Calculating evaluation metrics for the optimized model
mse_best = mean_squared_error(y_test, y_pred_best)
r2_best = r2_score(y_test, y_pred_best)

# Displaying the improvements in MSE and R^2
mse_improvement = mse - mse_best
r2_improvement = r2_best - r2

# Displaying the results
mse_best, r2_best, mse_improvement, r2_improvement #(0.2357857371499288, 0.8438206360172018, 1.404752938976063, 0.8017992644496518)


# Predicting the kp-index for the previously chosen Federicksburg k-index vector using the optimized model
sample_fredericksburg_scaled = scaler_X.transform(sample_fredericksburg_random)
predicted_kp_index_scaled_optimized = best_svr.predict(sample_fredericksburg_scaled)

# Transforming the predicted data back to the original scale
predicted_kp_index_optimized = scaler_y.inverse_transform(predicted_kp_index_scaled_optimized)

# Displaying the optimized prediction and actual kp-index
predicted_kp_index_optimized, actual_kp_index_random 
#(array([[3.31721741, 2.61251221,1.63346665, 2.87117847,1.03796742,
#0.98798582, 1.65333422, 0.67988594]]),
#array([[3.333, 2.667, 1.667, 2.667, 1. , 1. , 1.667, 0.667]]))



# Extracting data for the last 30 days
date_last_30 = date_columns['date'].tail(30)
X_last_30 = X.tail(30)
y_last_30 = y.tail(30)

# Predicting KP-index for the last 30 days
y_pred_last_30_scaled = best_svr.predict(scaler_X.transform(X_last_30))
y_pred_last_30 = scaler_y.inverse_transform(y_pred_last_30_scaled)


# Visualizing the actual vs. predicted KP-index for the last 30 days
plt.figure(figsize=(20, 15))

 # Looping through each KP-index column
for i in range(8):
    plt.subplot(4, 2, i+1)
    #Plotting the actual values
    plt.plot(date_last_30, y_last_30.iloc[:, i].values, 'r', label='Actual')
    #Plotting the predicted values
    plt.plot(date_last_30, y_pred_last_30[:, i], 'b--', label='Predicted')
    plt.title(f'KP-Index Column {i+1}')
    plt.xticks(rotation=45)
    plt.tight_layout() 
    if i == 0:
       plt.legend()

plt.suptitle('Actual vs Predicted KP-Index for Last 30 Days with Optimized Model', fontsize=16, y=1.05)
plt.show()


# Initializing the dataframe with the 'date' column
df_compact = pd.DataFrame()
df_compact['date'] = date_last_30.values

# Storing Fredericksburg values, predicted KP-index values, and actual KP-index values as lists
df_compact['federickburg'] = X_last_30.values.tolist()
df_compact['predicted_kp_index'] = y_pred_last_30.tolist()
df_compact['actual_kp_index'] = y_last_30.values.tolist()

# Rounding the values in 'predicted_kp_index' and 'actual_kp_index' columns to 2 decimal places
df_compact['predicted_kp_index'] = df_compact['predicted_kp_index'].apply(lambda x: [round(i, 2) for i in x])
df_compact['actual_kp_index'] = df_compact['actual_kp_index'].apply(lambda x: [round(i, 2) for i in x])

# Exporting the df_compact dataframe to a CSV file
csv_path = "pred_vs_act_last_30_days.csv"
df_compact.to_csv(csv_path, index=False)
df_compact.head(10) # Displaying the first few rows of the compact dataframe

# Calculating the evaluation metrics for the last 30 days' predictions
mae = mean_absolute_error(y_last_30, y_pred_last_30)
mse = mean_squared_error(y_last_30, y_pred_last_30)
r2 = r2_score(y_last_30, y_pred_last_30)
 
# Displaying the metrics
mae, mse, r2 #(0.06567479284845214, 0.010493205685711132, 0.9934506422188698)

My output images are:

enter image description here

enter image description here My question:

Is it possible to improve conciseness and complexity of my code? Or is there any alternative approach for this problem? Although I get the accuracy is almost 100% using the above algorithm.

Improve this question
\$\endgroup\$
2
  • \$\begingroup\$ You've made 10 revisions in a short amount of time. Do you consider the question finished now or are more incoming? \$\endgroup\$
    – Mast ♦
    Sep 20 at 20:06
  • \$\begingroup\$ @Mast I don't change one word till now except spelling mistaken words. \$\endgroup\$
    – S. M.
    Sep 20 at 20:12

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.